{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" }, "colab": { "name": "405_Adams Predictor Corrector Example.ipynb", "provenance": [], "include_colab_link": true } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "view-in-github", "colab_type": "text" }, "source": [ "\"Open" ] }, { "cell_type": "markdown", "metadata": { "id": "eHqI_Q92PMxU" }, "source": [ "# Adams Predictor Corrector\n", "The Adams Predictor corrector method is used to side step the issue of the Adams-Moulton method being implicit. It can also be used to estimate error when the solution is unknown." ] }, { "cell_type": "code", "metadata": { "id": "RiCsLfUZPMxV", "outputId": "93fc6c41-96b3-4794-899c-e6f39c9a7339", "colab": { "base_uri": "https://localhost:8080/", "height": 336 } }, "source": [ "from IPython.display import HTML\n", "HTML('')" ], "execution_count": 1, "outputs": [ { "output_type": "execute_result", "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "execution_count": 1 } ] }, { "cell_type": "markdown", "metadata": { "id": "fyV5F_mcPMxc" }, "source": [ "## Intial Value Problem\n", "\n", "The differential equation\n", "\\begin{equation}y^{'}=t-y, \\ \\ (0 \\leq t \\leq 2) \\end{equation}\n", "with the initial condition\n", "\\begin{equation}y(0)=1\\end{equation}\n", "will be used to illustrate the method." ] }, { "cell_type": "markdown", "metadata": { "id": "BgzGCYvjPMxd" }, "source": [ "## Python Libraries" ] }, { "cell_type": "code", "metadata": { "id": "tEQhmj5TPMxe" }, "source": [ "import numpy as np\n", "import math \n", "import pandas as pd\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt # side-stepping mpl backend\n", "import matplotlib.gridspec as gridspec # subplots\n", "import warnings\n", "\n", "warnings.filterwarnings(\"ignore\")\n" ], "execution_count": 2, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "yqb0HnEEPMxh" }, "source": [ "### Defining the function\n", "\\begin{equation} f(t,y)=t-y.$$" ] }, { "cell_type": "code", "metadata": { "id": "HZxaZLmAPMxi" }, "source": [ "def myfun_ty(t,y):\n", " return t-y" ], "execution_count": 3, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "EG5z-Dj-PMxk" }, "source": [ "## Discrete Interval\n", "Defining the step size $h$ from the interval range $a\\leq t \\leq b$ and number of steps $N$ \n", "\\begin{equation}h=\\frac{b−a}{N}.\\end{equation}\n", " \n", "This gives the discrete time steps,\n", "\\begin{equation}t_i=t_0+ih,\\end{equation}\n", "where $t0=a.$\n", "\n", "Here the interval is $0≤t≤2$ and number of step 4 \n", "\\begin{equation}h=\\frac{2−0}{4}=0.5.\\end{equation}\n", " \n", "This gives the discrete time steps,\n", "\\begin{equation}t_i=0+i0.5,\\end{equation}\n", "for $i=0,1,⋯,4.$" ] }, { "cell_type": "code", "metadata": { "id": "HNKrGexoPMxl", "outputId": "2bcd5f96-2c2c-4844-dae0-f2d59f800b4a", "colab": { "base_uri": "https://localhost:8080/", "height": 0 } }, "source": [ "# Start and end of interval\n", "b=2\n", "a=0\n", "# Step size\n", "N=4\n", "h=(b-a)/(N)\n", "t=np.arange(a,b+h,h)\n", "fig = plt.figure(figsize=(10,4))\n", "plt.plot(t,0*t,'o:',color='red')\n", "plt.xlim((0,2))\n", "plt.title('Illustration of discrete time points for h=%s'%(h))" ], "execution_count": 4, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "Text(0.5, 1.0, 'Illustration of discrete time points for h=0.5')" ] }, "metadata": {}, "execution_count": 4 }, { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmwAAAEICAYAAADiGKj0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAedklEQVR4nO3de7xldV3/8ddbRjHEuItcZ1DoAmVeTqBlNipysWQw6fcbQ0XUJgos62eFUl5ACs0C/akZKklAgJHmUCEXEfr9QIkzRAoqMozogCCXGS4jCqKf/ljrOHvO7HPOPu4zc9Y5vJ6Px3nM+q71XWt913d/997vsy5zUlVIkiSpux432w2QJEnS5AxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjapleS1Sf5/T7mS7D2bbZpIkg8n+fNZ2O/vJvl2knVJdhig/q1JDmyn35rko5u+lZtOkiOTXDLb7ZhMkj3b12eLzbCvlydZ3e7vWTOwvSuSvGEm2ibNNwY2Pab1BopNtP2PJ3nXkNvYIEgCVNUxVXXScK2bdjseD/wNcFBVbV1V905n/ar6i6rabF/GSRa1oXvBTK1fVedU1UEz18qZV1XfbF+fH0xVd9g+At4LHNfu779+zG3MuCTPTLIiyUPtv8+cpO4VSb7Xhs51SW7anG2VBmVgk2bREF+Us2Fn4InAjbPdEJhzfTdfLeTHHA+b6gxgkicAnwbOBrYDzgQ+3c6fyFjo3LqqfnpTtEsaloFNGsD4SzW9Z73SODXJXUkeSPKlJD+XZBlwJPAn7W/uF7b1b03yp0m+CHwnyYIkxye5JcmDSb6c5OVt3Z8FPgw8r93Gfe38Dc7cJfntJCuTrEmyPMmuPcsqyTFJbk5yX5IPJskEx7llktOSfKv9Oa2d91PA2JmH+5JcPsH6r07yjST3Jjlh3LJ3JDm7nX5ikrPbevcluTbJzu2y7ZP8fbv/tUn+pZ2/OMltbd/dCfx9ksf19N29ST6RZPt2l//R0951SZ7Xbud1Sb7SbvviJAsneNk3Wn/82c62b3+v7dsHk5yU5OlJrm7Hwid6g0KSX09yfXvMVyd5xgT7Htv27ydZleSeJH+V5HHtsscl+bO2r+9K8g9JtmmXbXDWrB27JyW5qm3jJUl2nOQY905yZZL72/2e36dtWyZZB2wB/HeSW9r5P9vu774kNyY5rGedjyf52yT/nuQ7wAsnOPSFE7R1UIuBBcBpVfVwVb0fCPCiaW5H6hQDmzS8g4AXAD8FbAP8L+DeqjodOAd4T/ub+8t61nkl8GvAtlX1KHAL8Cvt+u8Ezk6yS1V9BTgG+Hy7jW3H7zzJi4C/bPe7C/AN4Lxx1X4d+EXgGW29gyc4lhOA5wLPBH4B2B/4s6r6GrBfW2fbqtroyy/JvsDfAq8GdgV2AHafYD9Htce6R1vvGOC77bKzgK3a/T0FOLVnvacC29Oc2VkGvBE4HPjVdp9rgQ+2dV/Q096tq+rzSZYAbwV+A9gJ+H/AuRO0caP1J6h3MPAcmn77E+B04FXtsf0czWtNmnu8zgB+pz3mvwOWJ9lygu0CvBwYAZ4NLAFe185/bfvzQuBpwNbABybZzm8BR9P05xOAN09yjCcBl9Ccndod+L/jN9YGoa3b4i9U1dPTXDK/sF33KTSvzTlJes9Y/RZwMvBkYIPL/AO0lTYITvRzfFttP+CLteHfXfwi68dvP3/ZhtOrkiyepJ40awxs0vC+T/MF9DNAquorVXXHFOu8v6pWV9V3Aarqn6rqW1X1w6o6H7iZJiwN4kjgjKq6rqoeBt5Cc0ZuUU+dU6rqvqr6JvA5mkA20bZOrKq7qupumvD46gHbcQTwr1X1H207/hz44QR1v08TWvauqh9U1YqqeiDJLsChwDFVtbaqvl9VV/as90Pg7W1g+C5N0Duhqm5r9/kO4IhMfLn0GOAv29foUeAvgGdOcpZtEO+pqgeq6kbgBuCSqlpVVfcDFwFjN+MvA/6uqq5pj/lM4GGaoDeRd1fVmvZ1O402/NG8Tn/T7mcdzWu+dJLj/vuq+lrbZ59g4tcfmtdmIbBrVX2vqiYKVuM9lyY4nlJVj1TV5cC/9rQZ4NNVdVU7zr833bZW1baT/JzSVtsauH/cNu+neY/286c0oXc3mrB9YZKnD3jM0mZjYJOG1H4xfYDmzM5dSU5P8pNTrLa6t5DkNT2Xyu6jOTMz6KWgXWnOqo21Zx1wL80X0Jg7e6YfovlSm3Jb7fSuE9Ttt+6PjquqvtO2o5+zgIuB89pLn+9pz9DsAaypqrUTrHf3uC/6hcCnevrtK8APaO6362ch8L6e+mtoLpftNkH9QXy7Z/q7fcpjfb0Q+D+9Z4Vojney/u0dJ72vRb/XaQETH/egrz80ZwkD/Gd7WfN1k9TttSuwuqp6Q/o32LBvVzO16bS1n3XA+PffTwIP9qvcBugH218CzgSuAl46zX1Km5yBTRrMd2gu0415au/Cqnp/VT0H2Jfm0ugfjy2aYHs/mt+e3fkIcBywQ3vZ8waaL83JtjHmWzRhYGx7T6I5e3X7FOtNuS1gz3beIO6gCSBj7diqbcdG2jNn76yqfYFforlk+xqaL/Ttk2x06Xds1XHl1cCh4860PLGqbu9Td6z+74yr/xNVdfUA+xrWauDkcfveqqomuiQLPf3Jhq9Fv9fpUTYMi4PY6Bir6s6q+u2q2pXm8u2HMth/b/MtYI+x++x62tU7Dofq06x/krPfz1vbajcCz0g2uE/zGQz+cESx/r0ndYaBTRrM9cBvJNmq/fJ6/diCJL+Y5ID2DNF3gO+x/lLgt2kut0zmSTRfEne32zua5gzbmG8Du2fip9zOBY5O818ZbElzme+aqrp1OgfYs60/S7JTe7P322iethvEBcCvJ3l+29YTmeAzJskLk/x8micFH6C5DPfD9lLyRTQhYbskj0/ygn7baH0YOHnskmbb7iXtsrtpXoenjav/liT7tfW3SfKbE2y73/rD+AhwTDtWkuRJSX4tyUSX6gD+uO2HPYA/AMYeADgX+MMkeyXZmuY1P7+9zDsdGx1jkt9MMnbv4VqasTnRpe1e19CcEfuT9nVbDLyMje+n/LH1PMnZ7+cv2mpX0Jxl/f324Yjj2vkbPSiTZNskB6d5CGZBkiNp7uv7zEy1WZopBjZpMKcCj9CEpzNpHiYY85M0X8ZraS4B3Qv8VbvsY8C+7SWwf+m34ar6MvDXwOfb7f88zWWZMZfTnB24M8k9fda/jOZ+sX+mOcv1dGDpj3WU8C5glOYm7S8B17XzptTew3Us8I9tO9YCt01Q/ak0Ae8BmsuYV9JcJoXmnrnvA18F7gLeNMlu3wcsBy5J8iDwBeCAtj0P0dzgflXb/8+tqk8B76a5FPsAzZnMQyc4no3Wn7ITJlFVo8Bv01w+XwuspHlwYDKfBlbQ/MLwbzTjCZqHF86iecrz6zS/JLzxx2hTv2P8ReCaNE+BLgf+oKpWDbCtR2gC2qHAPcCHgNdU1Ven265htO04nOaM7X00D2oc3s4f+w+cL2qrP55mfN/dtvmNbd2vbc42S4PIhg/SSJK6IEkB+1TVytlui6TZ5xk2SZKkjjOwSZIkdZyXRCVJkjrOM2ySJEkdNyf/ePKOO+5YixYtmu1mSJIkTWnFihX3VNVOw2xjTga2RYsWMTo6OtvNkCRJmlKSb0xda3JeEpUkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeq4GQlsSQ5JclOSlUmO77N8yyTnt8uvSbJo3PI9k6xL8uaZaI8kSdJ8MnRgS7IF8EHgUGBf4JVJ9h1X7fXA2qraGzgVePe45X8DXDRsWyRJkuajmTjDtj+wsqpWVdUjwHnAknF1lgBnttMXAC9OEoAkhwNfB26cgbZIkiTNOzMR2HYDVveUb2vn9a1TVY8C9wM7JNka+FPgnVPtJMmyJKNJRu++++4ZaLYkSdLcMNsPHbwDOLWq1k1VsapOr6qRqhrZaaedNn3LJEmSOmLBDGzjdmCPnvLu7bx+dW5LsgDYBrgXOAA4Isl7gG2BHyb5XlV9YAbaJUmSNC/MRGC7FtgnyV40wWwp8Fvj6iwHjgI+DxwBXF5VBfzKWIUk7wDWGdYkSZI2NHRgq6pHkxwHXAxsAZxRVTcmOREYrarlwMeAs5KsBNbQhDpJkiQNIM2JrrllZGSkRkdHZ7sZkiRJU0qyoqpGhtnGbD90IEmSpCkY2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp4wxskiRJHWdgkyRJ6jgDmyRJUscZ2CRJkjrOwCZJktRxBjZJkqSOM7BJkiR1nIFNkiSp42YksCU5JMlNSVYmOb7P8i2TnN8uvybJonb+S5KsSPKl9t8XzUR7JEmS5pOhA1uSLYAPAocC+wKvTLLvuGqvB9ZW1d7AqcC72/n3AC+rqp8HjgLOGrY9kiRJ881MnGHbH1hZVauq6hHgPGDJuDpLgDPb6QuAFydJVf1XVX2rnX8j8BNJtpyBNkmSJM0bMxHYdgNW95Rva+f1rVNVjwL3AzuMq/MK4LqqengG2iRJkjRvLJjtBgAk2Y/mMulBk9RZBiwD2HPPPTdTyyRJkmbfTJxhux3Yo6e8ezuvb50kC4BtgHvb8u7Ap4DXVNUtE+2kqk6vqpGqGtlpp51moNmSJElzw0wEtmuBfZLsleQJwFJg+bg6y2keKgA4Ari8qirJtsC/AcdX1VUz0BZJkqR5Z+jA1t6TdhxwMfAV4BNVdWOSE5Mc1lb7GLBDkpXAHwFj//XHccDewNuSXN/+PGXYNkmSJM0nqarZbsO0jYyM1Ojo6Gw3Q5IkaUpJVlTVyDDb8C8dSJIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcTMS2JIckuSmJCuTHN9n+ZZJzm+XX5NkUc+yt7Tzb0py8EA7XLECFi2Cc86ZieZrPjrnnGaMPO5xjhVNzfGiQTlWNB3teHkOPGfYTS0YdgNJtgA+CLwEuA24NsnyqvpyT7XXA2urau8kS4F3A/87yb7AUmA/YFfgsiQ/VVU/mHLH3/gGLFvWTB955LCHofnknHOasfHQQ03ZsaLJOF40KMeKpmP8eBnSTJxh2x9YWVWrquoR4Dxgybg6S4Az2+kLgBcnSTv/vKp6uKq+DqxstzeYhx6Ct74VFi+Gs89eP2/xYjj//KZ8//1N+ZOfbMr33NOUL7ywKd95Z1P+zGea8urVTfmyy5ryqlVN+corm/JNNzXlq69uyjfc0JSvvbYpX399U77++qZ87bVN+YYbmvLVVzflm25qylde2ZRXrWrKl13WlFevbsqf+UxTvvPOpnzhhU35nnua8ic/2ZTvv78pn39+Ux4bIGef3ZS///2m/PGPN+UxH/kIHHjg+vKHPgSHHrq+/L73wWGHrS+/973wilesL59yCixdur580knwqletL7/tbXD00evLb3nL+g84gDe/GY49dn35TW9qfsYce2xTZ8yyZc02xhx9dLOPMa96VbPO+DfIQw/BCSc0bX/ve9fPP+yw5hjHHHpo0wdjDjyw6aMxixc3fQhNnzr25v7YO+GE/uPlmGOmP/ZOOml9eenSpo1jHHtzf+xNNFaOPXb2P/cce025S2Ov33gZwkwEtt2A1T3l29p5fetU1aPA/cAOA64LQJJlSUaTjG6wYPXqftX1WDb2Rhrvm9/cvO3Q3DDRuFi3bvO2Q9030ViZ6DNHj20z/J2TqhpuA8kRwCFV9Ya2/GrggKo6rqfODW2d29ryLcABwDuAL1TV2e38jwEXVdUFk+1zJKkfpbaFC+HWW4c6Bs0zixY1lyrGc6yoH8eLBuVY0XT0jJcRYLQqw2xuJs6w3Q7s0VPevZ3Xt06SBcA2wL0DrjuxrbaCk0+efos1v518cjM2ejlWNBHHiwblWNF09BsvQ5iJwHYtsE+SvZI8geYhguXj6iwHjmqnjwAur+bU3nJgafsU6V7APsB/DrTXhQvh9NO90VMbO/LIZmwsXAiJY0WTc7xoUI4VTUfveJkBQ18SBUjyUuA0YAvgjKo6OcmJwGhVLU/yROAs4FnAGmBpVa1q1z0BeB3wKPCmqrpoqv2NjIzU6OjoVNUkSZJmXZIVVTUy1DZmIrBtbgY2SZI0V8xEYPMvHUiSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOm6owJZk+ySXJrm5/Xe7Ceod1da5OclR7bytkvxbkq8muTHJKcO0RZIkab4a9gzb8cBnq2of4LNteQNJtgfeDhwA7A+8vSfYvbeqfgZ4FvDLSQ4dsj2SJEnzzrCBbQlwZjt9JnB4nzoHA5dW1ZqqWgtcChxSVQ9V1ecAquoR4Dpg9yHbI0mSNO8MG9h2rqo72uk7gZ371NkNWN1Tvq2d9yNJtgVeRnOWTpIkST0WTFUhyWXAU/ssOqG3UFWVpKbbgCQLgHOB91fVqknqLQOWAey5557T3Y0kSdKcNWVgq6oDJ1qW5NtJdqmqO5LsAtzVp9rtwOKe8u7AFT3l04Gbq+q0KdpxeluXkZGRaQdDSZKkuWrYS6LLgaPa6aOAT/epczFwUJLt2ocNDmrnkeRdwDbAm4ZshyRJ0rw1bGA7BXhJkpuBA9sySUaSfBSgqtYAJwHXtj8nVtWaJLvTXFbdF7guyfVJ3jBkeyRJkuadVM29q4sjIyM1Ojo6282QJEmaUpIVVTUyzDb8SweSJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHWcgU2SJKnjDGySJEkdZ2CTJEnqOAObJElSxxnYJEmSOs7AJkmS1HEGNkmSpI4zsEmSJHXcUIEtyfZJLk1yc/vvdhPUO6qtc3OSo/osX57khmHaIkmSNF8Ne4bteOCzVbUP8Nm2vIEk2wNvBw4A9gfe3hvskvwGsG7IdkiSJM1bwwa2JcCZ7fSZwOF96hwMXFpVa6pqLXApcAhAkq2BPwLeNWQ7JEmS5q1hA9vOVXVHO30nsHOfOrsBq3vKt7XzAE4C/hp4aKodJVmWZDTJ6N133z1EkyVJkuaWBVNVSHIZ8NQ+i07oLVRVJalBd5zkmcDTq+oPkyyaqn5VnQ6cDjAyMjLwfiRJkua6KQNbVR040bIk306yS1XdkWQX4K4+1W4HFveUdweuAJ4HjCS5tW3HU5JcUVWLkSRJ0o8Me0l0OTD21OdRwKf71LkYOCjJdu3DBgcBF1fV31bVrlW1CHg+8DXDmiRJ0saGDWynAC9JcjNwYFsmyUiSjwJU1Rqae9WubX9ObOdJkiRpAKmae7eDjYyM1Ojo6Gw3Q5IkaUpJVlTVyDDb8C8dSJIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeo4A5skSVLHGdgkSZI6zsAmSZLUcQY2SZKkjktVzXYbpi3Jg8BNs92OjtkRuGe2G9FB9kt/9kt/9svG7JP+7Jf+7Jf+frqqnjzMBhbMVEs2s5uqamS2G9ElSUbtk43ZL/3ZL/3ZLxuzT/qzX/qzX/pLMjrsNrwkKkmS1HEGNkmSpI6bq4Ht9NluQAfZJ/3ZL/3ZL/3ZLxuzT/qzX/qzX/obul/m5EMHkiRJjyVz9QybJEnSY4aBTZIkqeM6FdiSHJLkpiQrkxzfZ/mWSc5vl1+TZFHPsre0829KcvDmbPemNkC//FGSLyf5YpLPJlnYs+wHSa5vf5Zv3pZvWgP0y2uT3N1z/G/oWXZUkpvbn6M2b8s3nQH65NSe/vhakvt6ls3nsXJGkruS3DDB8iR5f9tvX0zy7J5l83WsTNUnR7Z98aUkVyf5hZ5lt7bzr5+J/66gSwbol8VJ7u95r7ytZ9mk77+5bIB++eOePrmh/TzZvl02L8dLkj2SfK79/r0xyR/0qTNzny1V1YkfYAvgFuBpwBOA/wb2HVfn94APt9NLgfPb6X3b+lsCe7Xb2WK2j2kz9ssLga3a6d8d65e2vG62j2EW++W1wAf6rLs9sKr9d7t2ervZPqbN0Sfj6r8ROGO+j5X22F4APBu4YYLlLwUuAgI8F7hmPo+VAfvkl8aOFTh0rE/a8q3AjrN9DLPUL4uBf+0zf1rvv7n2M1W/jKv7MuDy+T5egF2AZ7fTTwa+1ud7aMY+W7p0hm1/YGVVraqqR4DzgCXj6iwBzmynLwBenCTt/POq6uGq+jqwst3efDBlv1TV56rqobb4BWD3zdzG2TDIeJnIwcClVbWmqtYClwKHbKJ2bk7T7ZNXAudulpbNsqr6D2DNJFWWAP9QjS8A2ybZhfk7Vqbsk6q6uj1meOx8rgwyViYyzGdS502zXx4Tny1VdUdVXddOPwh8BdhtXLUZ+2zpUmDbDVjdU76NjQ/8R3Wq6lHgfmCHAdedq6Z7bK+nSfNjnphkNMkXkhy+KRo4Swbtl1e0p6EvSLLHNNedawY+rvay+V7A5T2z5+tYGcREfTdfx8p0jf9cKeCSJCuSLJulNs2m5yX57yQXJdmvnedYAZJsRRM8/rln9rwfL2lu0XoWcM24RTP22TJX/zSV+kjyKmAE+NWe2Qur6vYkTwMuT/Klqrpldlq42V0InFtVDyf5HZqzsy+a5TZ1xVLggqr6Qc+8x/JY0QSSvJAmsD2/Z/bz27HyFODSJF9tz8A8FlxH815Zl+SlwL8A+8xym7rkZcBVVdV7Nm5ej5ckW9ME1DdV1QObaj9dOsN2O7BHT3n3dl7fOkkWANsA9w647lw10LElORA4ATisqh4em19Vt7f/rgKuoPkNYD6Ysl+q6t6evvgo8JxB152jpnNcSxl3yWIej5VBTNR383WsDCTJM2jeO0uq6t6x+T1j5S7gU8yfW1CmVFUPVNW6dvrfgccn2ZHH+FjpMdlny7wbL0keTxPWzqmqT/apMnOfLbN9017PjXkLaG6624v1N2zuN67OsWz40MEn2un92PChg1XMn4cOBumXZ9Hc7LrPuPnbAVu20zsCNzNPboIdsF926Zl+OfCFdnp74Ott/2zXTm8/28e0OfqkrfczNDcB57EwVnqOcRET30j+a2x4Y/B/zuexMmCf7ElzP/AvjZv/JODJPdNXA4fM9rFsxn556th7hyZ4fLMdNwO9/+byz2T90i7fhuY+tyc9FsZL+7r/A3DaJHVm7LOlM5dEq+rRJMcBF9M8bXNGVd2Y5ERgtKqWAx8DzkqykmZQLG3XvTHJJ4AvA48Cx9aGl3rmrAH75a+ArYF/ap7B4JtVdRjws8DfJfkhzdnUU6rqy7NyIDNswH75/SSH0YyJNTRPjVJVa5KcBFzbbu7E2vD0/Zw0YJ9A8745r9pPjda8HSsASc6lebpvxyS3AW8HHg9QVR8G/p3maa6VwEPA0e2yeTlWYKA+eRvNPcIfaj9XHq2qEWBn4FPtvAXAP1bVZzb7AWwiA/TLEcDvJnkU+C6wtH0v9X3/zcIhbBID9As0vxhfUlXf6Vl1Po+XXwZeDXwpyfXtvLfS/LIz458t/mkqSZKkjuvSPWySJEnqw8AmSZLUcQY2SZKkjjOwSZIkdZyBTZIkqeMMbJIkSR1nYJMkSeq4/wGMoExHZjwsOQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "f0oBK4BePMxn" }, "source": [ "## Exact Solution\n", "The initial value problem has the exact solution\n", "\\begin{equation}y(t)=2e^{-t}+t-1.\\end{equation}\n", "The figure below plots the exact solution." ] }, { "cell_type": "code", "metadata": { "id": "2fvucY9iPMxn", "outputId": "e336ecd6-f139-4bfd-f78d-0d5a490c3725", "colab": { "base_uri": "https://localhost:8080/", "height": 312 } }, "source": [ "IC=1 # Intial condtion\n", "y=(IC+1)*np.exp(-t)+t-1\n", "fig = plt.figure(figsize=(6,4))\n", "plt.plot(t,y,'o-',color='black')\n", "plt.title('Exact Solution ')\n", "plt.xlabel('time')" ], "execution_count": 5, "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "Text(0.5, 0, 'time')" ] }, "metadata": {}, "execution_count": 5 }, { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEWCAYAAAB2X2wCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhU9dnG8e8DsguCLBajJGJ5iyiKktpKRcAVkIIbAgaEiCxWVKTghoqVRpGyK4tBwYoREnGFIpRdXyhIBMFdURbDIggiFCQQ8nv/mInvGAJJYGbOLPfnuuZi5pwzc+6cDM/88pwz55hzDhERiX5lvA4gIiLBoYIuIhIjVNBFRGKECrqISIxQQRcRiREq6CIiMUIFXeQEmNlGM7v6BJ9bz8z+a2Zlg51L4psKuoSFvwD+7C9kBbfnQri+lmaWU8wyZ5nZ62b2g5n9ZGafmFmPEGT5VfF3zm12zp3qnDsS7HVJfDvF6wASV/7snFvgdYgA04C1QCKQCzQGfuNpIpGToBG6eM7MJprZ6wGPnzGzheZTw8xmm9lOM/vRf/+sgGVPN7OpZrbVP/8tM6sCvAucGfDXwJlFrPr3wEvOuf3OuTzn3Brn3LsBr93ezD41sz1mtsTMzjtG/pfM7O8Bj3/568DMpgH1gFn+HA+YWZKZOTM7xb/MmWb2jpntNrP1ZtYr4LWeMLMsM3vZzPb58ySf6LaW2KaCLpHgr0BjM+thZs2BnkB35zsvRRlgKr5RdD3gZyCwVTMNqAycD9QBRjvn9gNtgK3+1sapzrmtRax3BTDezDqbWb3AGWb2P8B0oD9QG5iDryiXL80P5pzrBmzG99fJqc654UUsNgPIAc4EbgGeMrMrA+a39y9THXin0M8v8gsVdAmnt/yj3YJbLwDn3AGgGzAKeAW4xzmX45+3yzn3unPugHNuH5AGtAAws7r4Cndf59yPzrnDzrmlpcjTEXgfeAzYYGYfmdnv/fM6Af9yzs13zh0GRgCVgGYnuQ1+xczOBv4EPOicO+ic+wh4Abg9YLH/dc7N8ffcpwEXBTODxA4VdAmnG5xz1QNukwtmOOdWAt8CBmQVTDezymb2vJltMrO9wHtAdf8RImcDu51zP55IGP+HwEPOufOBM4CP8H3oGL7R8qaAZfOB74CEE1nXcZyJ72fYFzBtU6H1bA+4fwCoWNCuEQmkgi4RwczuBioAW4EHAmb9Ffgd8AfnXDXgioKn4Cuwp5tZ9SJeslSnEXXO/YBvFH4mcLo/R2JAPsP3AbKliKfvx9f2KVB4x+rxsmzF9zNUDZhW7xjrETkuFXTxnL9f/XegK77WywNm1sQ/uyq+vvkeMzsdGFLwPOfcNnw7Pyf4d56WM7OCgv89UNPMTjvOep8xswvM7BR/Qb0LWO+c24Xvr4TrzewqMyuH74MlF1hexEt9BLT176D9Db6+e6DvgfpFZXDOfed/zafNrKKZXYhvH8Irx8otciwq6BJOBUd6FNze9LcOXgGecc6tdc59DTwCTDOzCsAYfL3rH/DtxJxb6DW7AYeBL4Ad+Iupc+4LfDs1v/X364s6yqUy8CawB1+7JxHfDkicc1/i+4B51r/uP+PbsXmoiNcpOPxxI/BvILPQ/KeBR/05Bhbx/C5AEr7R+pvAkAg7vFOihOkCFyIisUEjdBGRGKGCLiISI1TQRURihAq6iEiM8OzLCbVq1XJJSUlerV5EJCp9+OGHPzjnahc1z7OCnpSURHZ2tlerFxGJSma26Vjz1HIREYkRKugiIjFCBV1EJEaooIuIxAgVdBGRGKGCLiISJhkZGSQlJVGmTBmSkpLIyMgI6uvrJPkiImGQkZFB7969OXDgAACbNm2id+/eAKSkpARlHRqhi4iEweDBg38p5gUOHDjA4MGDg7YOFXQRkTDYvHlzqaafCBV0EZEwSEgo+nK09erVC9o6VNBFRMLg3HPPPWpa5cqVSUtLC9o6VNBFREJs7ty5LF26lPbt25OYmIiZkZiYSHp6etB2iIKOchERCam9e/fSq1cvzjvvPDIzM6lYsWLI1qWCLiISQoMGDWLr1q0sX748pMUc1HIREQmZBQsWkJ6ezl//+lf+8Ic/hHx95pwL+UqKkpyc7HQ+dBGJVfv27aNx48ZUqFCBjz76iEqVKgXldc3sQ+dcclHz1HIREQmBBx98kM2bN/O///u/QSvmxVHLRUQkyBYvXszEiRPp378/zZo1C9t61XIREQmi/fv307hxY8qWLcvatWupXLlyUF9fLRcRkTB5+OGH2bhxI0uXLg16MS+OWi4iIkHy/vvv8+yzz9KvXz+aN28e9vWr5SIiEgQHDhzgoosu4siRI3z88cdUqVIlJOtRy0VEJMQeffRR1q9fz6JFi0JWzIujlouIyElavnw5Y8aM4a677qJVq1ae5VBBFxE5CT///DOpqanUq1ePZ555xtMsarmIiJyEIUOG8NVXXzF//nyqVq3qaRaN0EVETtDKlSsZOXIkvXr14uqrr/Y6jgq6iMiJOHjwIKmpqSQkJDBixAiv4wBquYiInJAnn3ySzz//nLlz51KtWjWv4wAaoYuIlFp2djbDhw/njjvu4LrrrvM6zi+KLehmNsXMdpjZJ8eYn2Jm68zsYzNbbmYXBT+miEhkyM3NJTU1lTPOOIORI0d6HedXSjJCfwlofZz5G4AWzrnGwFAgPQi5REQiUlpaGp988gnp6elUr17d6zi/UmwP3Tn3npklHWf+8oCHK4CzTj6WiEjkWbNmDU899RS33347119/vddxjhLsHnpP4N1jzTSz3maWbWbZO3fuDPKqRURC59ChQ6SmplK7dm1Gjx7tdZwiBe0oFzNrha+gX36sZZxz6fhbMsnJyd6cFUxE5AQMGzaMtWvX8tZbb3H66ad7HadIQSnoZnYh8ALQxjm3KxivKSISKdatW8fQoUO57bbb6NChg9dxjumkWy5mVg94A+jmnPvq5COJiESOw4cPk5qayumnn864ceO8jnNcxY7QzWw60BKoZWY5wBCgHIBzbhLwOFATmGBmAHnHOleviEi0GT58OKtXr+b111+nZs2aXsc5Ll3gQkTkGD799FMuueQSbrjhBjIzM72OAxz/Ahf6pqiISBHy8vJITU2lWrVqPPfcc17HKRGdy0VEpAgjR45k1apVZGZmUrt2ba/jlIhG6CIihXzxxRcMGTKEm266iY4dO3odp8RU0EVEAhw5coTU1FSqVKnChAkT8B/sERXUchERCTBmzBhWrFhBRkYGZ5xxhtdxSkUjdBERv6+++opHH32U9u3b06VLF6/jlJoKuogIvlbLHXfcQcWKFZk0aVJUtVoKqOUiIgI899xzLFu2jH/+85/UrVvX6zgnRCN0EYl769ev5+GHH6Zt27Z069bN6zgnTAVdROJafn4+PXv2pHz58qSnp0dlq6WAWi4iEtcmTpzIe++9x4svvkhCQoLXcU6KRugiErc2bNjAgw8+yHXXXUdqaqrXcU6aCrqIxCXnHHfeeSdlypSJ+lZLAbVcRCQupaens2jRIp5//nnq1avndZyg0AhdROLOpk2bGDhwIFdffTW9evXyOk7QqKCLSFxxztGrVy+cc0yePDkmWi0F1HIRkbgyZcoU5s+fz/jx40lKSvI6TlBphC4icSMnJ4cBAwbQsmVL+vbt63WcoFNBF5G44Jyjd+/e5OXl8eKLL1KmTOyVP7VcRCQuvPzyy7z77ruMHTuW+vXrex0nJGLvI0pEpJCtW7fSv39/mjdvTr9+/byOEzIq6CIS05xz9OnTh4MHD8Zsq6WAWi4iEtMyMjKYPXs2o0aNokGDBl7HCanY/agSkbi3fft27r33Xi677DLuvfder+OEnAq6iMQk5xx33XUXBw4cYMqUKZQtW9brSCGnlouIxKTMzEzeeusthg8fTsOGDb2OExYaoYtIzNmxYwf9+vXj0ksvZcCAAV7HCZtiC7qZTTGzHWb2yTHmNzSz/5hZrpkNDH5EEZHSufvuu9m3bx9Tp06Ni1ZLgZKM0F8CWh9n/m7gXmBEMAKJiJyMmTNnMnPmTJ544gkaNWrkdZywKragO+few1e0jzV/h3NuFXA4mMFERErrhx9+4C9/+QtNmzZl0KBBXscJu7DuFDWz3kBvIGZOKC8ikeOee+5hz549LFy4kFNOib9jPsK6U9Q5l+6cS3bOJdeuXTucqxaRGPfWW28xY8YMHnvsMRo3bux1HE/oKBcRiXq7d++mb9++NGnShIceesjrOJ6Jv79JRCTm3HfffezatYu5c+dSrlw5r+N4ptiCbmbTgZZALTPLAYYA5QCcc5PM7DdANlANyDez/kAj59zekKUWEfGbNWsWr7zyCo8//jhNmjTxOo6nzDnnyYqTk5Nddna2J+sWkdjw448/cv7551OrVi2ys7MpX76815FCzsw+dM4lFzVPLRcRiVoDBgxgx44dzJ49Oy6KeXG0U1REotK7777LSy+9xIMPPsgll1zidZyIoJaLiESdn376ifPPP5/TTjuN1atXU6FCBa8jhY1aLiISUwYOHMi2bdt4880346qYF0ctFxGJKv/+97954YUXGDhwIL///e+9jhNR1HIRkaixb98+LrjgAipXrsyaNWuoWLGi15HCTi0XEYkJDzzwAN999x3Lli2Ly2JeHLVcRCQqLFq0iEmTJjFgwAAuu+wyr+NEJLVcRCTi/fe//6Vx48aUK1eOtWvXUqlSJa8jeUYtFxGJag8//DCbNm3ivffei+tiXhy1XEQkoi1dupTnnnuOe++9l8svv9zrOBFNLRcRiVj79+/noosuwjnHunXrqFKliteRPKeWi4hEpcGDB/PNN9+wePFiFfMSUMtFRCLSsmXLGDduHHfffTctW7b0Ok5UUEEXkYjz888/k5qaSmJiIsOGDfM6TtRQy0VEIs5jjz3G119/zcKFCzn11FO9jhM1NEIXkYiyYsUKRo8eTZ8+fbjyyiu9jhNVVNBFJGIcPHiQ1NRUEhISGD58uNdxoo5aLiISMZ544gm++OIL5s2bR7Vq1byOE3U0QheRiLBq1Sr+8Y9/0LNnT6699lqv40QlFXQR8Vxubi49evSgbt26jBw50us4UUstFxHx3NChQ/nss8/417/+xWmnneZ1nKilEbqIeGr16tUMGzaM7t2707ZtW6/jRDUVdBHxzKFDh+jRowd16tRh9OjRXseJemq5iIhnnnrqKT7++GPeeecdatSo4XWcqKcRuoh4Yu3ataSlpZGSksKf//xnr+PEhGILuplNMbMdZvbJMeabmY0zs/Vmts7MLgl+TJ+MjAySkpIoU6YMSUlJZGRkhGpVIhJChw8fpkePHtSsWZOxY8d6HSdmlGSE/hLQ+jjz2wAN/LfewMSTj3W0jIwMevfuzaZNm3DOsWnTJnr37q2iLhKFnnnmGT766CMmTpxIzZo1vY4TM0p0gQszSwJmO+cuKGLe88AS59x0/+MvgZbOuW3He83SXuAiKSmJTZs2HTU9MTGRjRs3lvh1RMRbH3/8MU2bNuXmm29m+vTpXseJOse7wEUweugJwHcBj3P804oK0tvMss0se+fOnaVayebNm0s1XUQiT15eHqmpqVSvXp1nn33W6zgxJ6w7RZ1z6c65ZOdccu3atUv13Hr16pVquohEnhEjRvDhhx8yfvx4atWq5XWcmBOMgr4FODvg8Vn+aUGVlpZG5cqVj5p+3333BXtVIhICn332GUOGDOGWW26hY8eOXseJScEo6O8At/uPdvkj8FNx/fMTkZKSQnp6OomJiZgZCQkJVKpUiVdeeYXc3Nxgr05EgujIkSPccccdVK1alfHjx3sdJ2aV5LDF6cB/gN+ZWY6Z9TSzvmbW17/IHOBbYD0wGfhLqMKmpKSwceNG8vPzycnJYfr06axevZqBAweGapUiEgSjR49m5cqVPPvss9SpU8frODGrREe5hEJpj3I5lgEDBjB69GhmzpzJzTffHIRkIhJMX375JRdddBFt2rThjTfewMy8jhTVjneUS9QX9EOHDtG8eXO++OIL1qxZQ/369YOQTkSC4ciRI1xxxRV8/vnnfPbZZ/zmN7/xOlLUC/Vhi54qX748M2bMoEyZMnTq1En9dJEIMm7cOJYvX864ceNUzMMg6gs6wDnnnMPUqVPJzs7mgQce8DqOiADr169n8ODBtGvXjpSUFK/jxIWYKOgAN9xwA/fddx/jxo3jzTff9DqOSFzLz8/njjvuoHz58kyaNEl98zCJmYIOMHz4cJKTk0lNTWXDhg1exxGJW+PHj+f9999nzJgxJCQU+cVxCYGYKujly5cnMzMTgE6dOnHo0CGPE4nEn2+//ZaHHnqINm3a0L17d6/jxJWYKugA9evXZ8qUKaxatYoHH3zQ6zgicSU/P5+ePXtStmxZnn/+ebVawizmCjrATTfdxD333MOYMWN4++23vY4jEjeef/55lixZwqhRozj77LOLf4IEVdQfh34subm5/OlPf+Kbb75hzZo1JCUlhWxdIgIbN26kcePGXHbZZcybN0+j8xCJ6ePQj6VChQpkZmaSn59P586d1U8XCSHnHL169QJg8uTJKuYeidmCDnDuuefy4osvsnLlSh5++GGv44jEnMDLQi5YsICOHTuSmJjoday4FdMFHeCWW27h7rvvZtSoUbzzzjtexxGJGYGXhSyQmZmpy0J6KGZ76IEOHjxIs2bN2LhxI2vWrNEIQiQIdFlIb8RlDz1QxYoVycrKIi8vj86dO3P48GGvI4lEPV0WMvLERUEH+O1vf8sLL7zAihUreOSRR7yOIxLV8vLyqFSpUpHzdFlI78RNQQe49dZbueuuuxgxYgSzZ8/2Oo5IVDp8+DBdunThwIEDlCtX7lfzKleuTFpamkfJJK4KOsCoUaNo0qQJ3bt357vvvvM6jkhUOXToEJ06dWLmzJmMGjWKqVOn/nJZyMTERNLT03VmRQ/FxU7Rwr7++msuueQSLrzwQpYsWXLUKENEjpabm0vHjh2ZNWsW48aN45577vE6UlyK+52ihTVo0IDJkyezfPlyHn30Ua/jiES8gwcPcuONNzJr1iwmTpyoYh6h4rKgA3Tu3Jk+ffowfPhw5syZ43UckYh14MAB2rdvz9y5c5k8eTJ9+/Yt/kniibgt6OC7EvmFF17I7bffTk5OjtdxRCLO/v37adeuHQsWLGDq1KnceeedXkeS44jrgl6pUiWysrLIzc2lc+fO5OXleR1JJGLs27ePtm3bsnTpUqZNm6Zzm0eBuC7oAL/73e94/vnnWbZsGY899pjXcUQiwt69e2ndujXLli3j1Vdf1ZErUSLuCzrAbbfdRq9evRg2bBhz5871Oo6Ip/bs2cO1117LBx98QGZmJp06dfI6kpSQCrrf2LFjady4Md26dWPLli1exxHxxO7du7nmmmtYvXo1M2fO5Oabb/Y6kpSCCrpfQT/9559/pkuXLuqnS9zZtWsXV111FevWreONN96gQ4cOXkeSUipRQTez1mb2pZmtN7OHipifaGYLzWydmS0xs7OCHzX0GjZsyKRJk3j//fcZMmSI13FEwmbnzp20atWKzz//nLfffpt27dp5HUlOQLEF3czKAuOBNkAjoIuZNSq02AjgZefchcCTwNPBDhouXbt2pWfPnjz99NPMmzfP6zgiIff999/TqlUr1q9fz+zZs2ndurXXkeQElWSEfimw3jn3rXPuEDADKPy3WCNgkf/+4iLmR5Vx48Zx/vnn061bN7Zu3ep1HJGQ2bZtGy1btmTDhg3MmTOHq6++2utIchJKUtATgMCzWOX4pwVaC9zkv38jUNXMap58PG9UrlyZrKws9u/fz2233aZ+usSknJwcWrRoQU5ODnPnzqVly5ZeR5KTFKydogOBFma2BmgBbAGOFF7IzHqbWbaZZe/cuTNIqw6N8847j4kTJ7J06VL+9re/eR1HJKg2b95MixYt2L59O/PmzaN58+ZeR5IgKElB3wKcHfD4LP+0XzjntjrnbnLOXQwM9k/bU/iFnHPpzrlk51xy7dq1TyJ2eNx+++2kpqaSlpbG/PnzvY4jEhQbN26kRYsW7Nq1iwULFtCsWTOvI0mQlKSgrwIamNk5ZlYe6Az86mrLZlbLzApe62FgSnBjeufZZ5/lvPPOIyUlhW3btnkdR+SkfPPNN1xxxRX89NNPLFy4kEsvvdTrSBJExRZ051we0A+YB3wOZDnnPjWzJ82svX+xlsCXZvYVcAYQM5csqVKlCq+99tov/fQjR47qJIlEha+++ooWLVpw4MABFi1aRNOmTb2OJEEWlxe4OBEvvfQSqampPP744+qpS9T54osvuPLKK8nLy2PhwoU0btzY60hygnSBiyDo0aMH3bt3Z+jQoSxcuNDrOCIl9umnn9KiRQvy8/NZsmSJinkMU0EvhfHjx9OwYUNSUlLYvn2713FEirVu3TpatmxJ2bJlWbJkCY0aFf5OoMQSFfRSqFKlCllZWezdu5eUlBT10yWirVmzhlatWlGxYkWWLl1Kw4YNvY4kIaaCXkoXXHABzz33HIsWLeLvf/+713FEipSdnc2VV17JqaeeytKlS2nQoIHXkSQMVNBPQGpqKt26deNvf/sbixYtKv4JImG0cuVKrr76amrUqMHSpUupX7++15EkTFTQT4CZMWHCBH73u9+RkpLC999/73UkEQCWLVvGNddcQ61atVi6dClJSUleR5IwUkE/QaeeeipZWVns2bOHrl27qp8unnvvvfe47rrrqFu3LkuXLuXss88u/kkSU1TQT0Ljxo159tlnWbBgAU899ZTXcSSOLVq0iDZt2lCvXj2WLFlCQkLh8+dJPFBBP0k9e/YkJSWFJ554giVLlngdR+LQv//9b66//nrq16/P4sWLqVu3rteRxCMq6CfJzJg4cSK//e1vue2229ixY4fXkSSOzJkzh/bt2/M///M/LFq0iDPOOMPrSOIhFfQgqFq1Kq+99ho//vgjXbt2JT8/3+tIEgdmzZrFjTfeSKNGjVi0aBHRcAZTCS0V9CC58MILGTt2LPPnz+fpp6P2CnwSJd58801uuukmLrroIhYuXEjNmlF7PRkJIhX0IOrVqxddunTh8ccf57333vM6jsSo1157jY4dO5KcnMz8+fOpUaOG15EkQqigB5GZ8fzzz3PuuefSpUsXIv2qTBJ9pk+fTpcuXfjjH//IvHnzOO2007yOJBFEBT3IqlatSlZWFrt27aJbt27qp0vQTJs2ja5du3L55Zczd+5cqlWr5nUkiTAq6CHQpEkTxowZw7x583jmmWe8jiMxYMqUKXTv3p1WrVoxZ84cTj31VK8jSQRSQQ+RPn360KlTJx577DHef/99r+NIFEtPT6dnz55ce+21zJo1i8qVK3sdSSKUCnqImBnp6emcc845dOnShR9++MHrSBKFxo8fT58+fbj++ut56623qFSpkteRJIKpoIdQtWrVyMrKYufOndx+++3qp0upjB07ln79+tGhQwdef/11Klas6HUkiXAq6CF28cUXM3r0aN59913+8Y9/eB1HosSIESPo378/N910E1lZWVSoUMHrSBIFVNDD4K677qJjx44MHjyYZcuWeR1HItzTTz/NoEGDuPXWW5kxYwbly5f3OpJECRX0MDAzJk+eTGJiIp07d2bXrl1eR5II9eSTT/LII49w2223kZGRQbly5byOJFFEBT1MTjvtNLKystixYwfdu3dXP11+xTnH448/zpAhQ+jevTsvv/wyp5xyitexJMqooIdR06ZNGTlyJP/6178YOXKk13EkQjjneOSRRxg6dCg9e/ZkypQplC1b1utYEoVU0MPs7rvv5uabb+bhhx9m+fLlXscRjznnGDRoEMOGDaNv376kp6dTpoz+W8qJ0TsnzMyMF198kXr16qmfHuecc9x///2MHDmSfv36MWHCBBVzOSl693igoJ++fft2evTogXPO60gSZvn5+fTr14+xY8dy//33M27cOMzM61gS5UpU0M2stZl9aWbrzeyhIubXM7PFZrbGzNaZWdvgR40tycnJjBgxgtmzZzNq1Civ40gY5efn07dvXyZMmMADDzzAyJEjVcwlOJxzx70BZYFvgPpAeWAt0KjQMunAXf77jYCNxb1u06ZNXbzLz893N954ozvllFPcf/7zH6/jSBjk5eW51NRUB7jBgwe7/Px8ryNJlAGy3THqaklG6JcC651z3zrnDgEzgA6FPxeAgnN5ngZsPYnPmLhhZkyZMoWzzjqLTp06sXv3bq8jSQgdOXKE1NRUpk6dyhNPPMHQoUM1MpegKklBTwC+C3ic458W6Amgq5nlAHOAe4p6ITPrbWbZZpatiz/4VK9enczMTLZt20Zqaqr66TEqLy+Prl27Mm3aNP7+978zZMgQFXMJumDtFO0CvOScOwtoC0wzs6Ne2zmX7pxLds4l64K2/+/SSy9l+PDhvPPOO4wZM8brOBJkhw8fpkuXLsyYMYNnnnmGwYMHex1JYlRJCvoW4OyAx2f5pwXqCWQBOOf+A1QEagUjYLy477776NChAw8++CAffPCB13EkSA4dOkSnTp2YOXMmo0aN4oEHHvA6ksSwkhT0VUADMzvHzMoDnYF3Ci2zGbgKwMzOw1fQ1VMpBTNj6tSpnHnmmdx66638+OOPXkeSk5Sbm8stt9zCm2++ybhx47j//vu9jiQxrtiC7pzLA/oB84DPgSzn3Kdm9qSZtfcv9legl5mtBaYDPZyawaVWo0YNMjMz2bJlC3fccYf66VHs4MGD3HjjjcyaNYuJEydyzz1F7lYSCa5jHf4S6psOWzy2UaNGOcCNGTPG6yhyAvbv3++uueYaZ2Zu8uTJXseRGMNJHrYoYda/f3/at2/PoEGDWLVqlddxpBT2799Pu3btWLBgAVOnTuXOO+/0OpLEERX0CFTQT69bty633nore/bs8TqSlMC+ffto27YtS5cuZdq0aXTv3t3rSBJnVNAj1Omnn05mZiY5OTn07NlT/fQIt3fvXlq3bs2yZct49dVXSUlJ8TqSxCEV9Aj2xz/+kWHDhvHGG2/w3HPPeR1HjmHPnj1ce+21fPDBB2RmZtKpUyevI0mcUkGPcAMGDPoA6UEAAAyGSURBVKBdu3YMHDiQ7Oxsr+NIIbt37+aaa65h9erVzJw5k5tvvtnrSBLHVNAjnJnx0ksvccYZZ9CpUyd++uknryOJ365du7jqqqtYt24db7zxBh06FD7FkUh4qaBHgZo1azJjxgw2bdrEnXfeqX56BNi5cyetWrXi888/5+2336Zdu3ZeRxJRQY8WzZo14+mnn2bmzJlMmDDB6zhxbfv27bRs2ZL169cze/ZsWrdu7XUkEUAFPar89a9/pW3btgwYMIDVq1d7HScubd26lZYtW7Jx40bmzJnD1Vdf7XUkkV+ooEeRMmXK8M9//pM6depw6623snfvXq8jxZWcnBxatmzJli1bmDt3Li1btvQ6ksivqKBHmVq1ajFjxgw2btxIr1691E8Pk02bNtGiRQu2b9/OvHnzaN68udeRRI6igh6F/vSnP5GWlkZWVhaTJk3yOk7M27BhAy1atGDXrl0sWLCAZs2aeR1JpEgq6FFq0KBBtGnThv79+7NmzRqv48Ss9evX06JFC/bu3cvChQu59NJLvY4kckwq6FGqoJ9eu3Zt9dND5KuvvqJFixYcOHCARYsW0bRpU68jiRyXCnoUq127NtOnT2fDhg307t1b/fQgyMjIICkpiTJlynDeeeexb98+Fi9eTJMmTbyOJlIsFfQo17x5c4YOHUpmZibp6elex4lqGRkZ9O7dm02bNuGcIz8/n8OHD7Nu3Tqvo4mUiHk1qktOTnY6N0lw5Ofn07ZtW5YsWcKKFSs0mizGwYMH2bJly1G39PR0Dhw4cNTyiYmJbNy4MfxBRYpgZh8655KLnKeCHht27NjBxRdfTJUqVfjwww+pWrWq15HCzjnHrl27finQOTk5RRbu3bt3H/XcKlWqsH///iJf18zIz88PdXyREjleQT8l3GEkNOrUqcP06dNp1aoVffr0ISMjAzPzOlbQ5ObmsnXr1iILdMFt69at5Obm/up5ZkadOnVISEggKSmJyy+/nISEhKNu1apV45xzzmHTpk1HrbtevXrh+jFFTooKegy54oorePLJJ3n00Udp1aoVvXr18jpSsZxz/Pjjj8cs0gWj7B9++OGo51aqVOmXgnzZZZcVWajr1q1LuXLlSpQlLS2N3r17/6rtUrlyZdLS0oL284qEklouMSY/P5/WrVvz/vvvs3LlSi688ELPshw+fJht27Yds/VRcDt48OBRz61du3aRBTohIYGzzjqLhIQEqlevHvS/QjIyMhg8eDCbN2+mXr16pKWl6epDElHUQ48zO3bsoEmTJlStWpXs7Oyg99Odc/z000/HLdJbtmxhx44dRx1KWaFCBc4888yjinPhUXWFChWCmlkkVqiHHmfq1KnDq6++ylVXXUXr1q3Jycnhu+++K9GIMy8vj+3btx+z9VFwK+pokJo1a/5SlC+55JIiR9c1a9aMqd6+SCTRCD2G3XLLLbz++uu/mlaxYkX69+9Pw4YNixxVf//990cd0VGuXLlfjaqLup155plUqlQpnD+eSFzSCD1OrVq16qhpBw8eZNiwYb88rl69+i9FuXHjxkW2QWrVqkWZMvoOmkikU0GPYd99912R082ML7/8koSEBCpXrhzmVCISKiUadplZazP70szWm9lDRcwfbWYf+W9fmdme4EeV0jrW8dP16tWjQYMGKuYiMabYgm5mZYHxQBugEdDFzBoFLuOcu98518Q51wR4FngjFGGldNLS0o4q2jquWiR2lWSEfimw3jn3rXPuEDAD6HCc5bsA04MRTk5OSkoK6enpJCYmYmYkJiaSnp6u46pFYlRJeugJQGAzNgf4Q1ELmlkicA6w6OSjSTCkpKSogIvEiWAfutAZmOmcO1LUTDPrbWbZZpa9c+fOIK9aRCS+laSgbwHODnh8ln9aUTpznHaLcy7dOZfsnEuuXbt2yVOKiEixSlLQVwENzOwcMyuPr2i/U3ghM2sI1AD+E9yIIiJSEsUWdOdcHtAPmAd8DmQ55z41syfNrH3Aop2BGU7XQRMR8USJvljknJsDzCk07fFCj58IXiwRESktz87lYmY7gaOvJlAytYCjT5DtvUjNBZGbTblKR7lKJxZzJTrnitwJ6VlBPxlmln2sk9N4KVJzQeRmU67SUa7SibdcOuOSiEiMUEEXEYkR0VrQ070OcAyRmgsiN5tylY5ylU5c5YrKHrqIiBwtWkfoIiJSiAq6iEiMiLiCXoKLaVQws0z//JVmlhQw72H/9C/N7Low5xpgZp+Z2TozW+g/82TBvCMBFwA56rQJIc7Vw8x2Bqz/zoB53c3sa/+te5hzHfOiKCHeXlPMbIeZfXKM+WZm4/y515nZJQHzQrm9isuV4s/zsZktN7OLAuZt9E//yMyCeqHeEuRqaWY/Bfy+Hg+Yd9z3QIhzDQrI9In/PXW6f15ItpeZnW1mi/114FMzu6+IZUL7/nLORcwNKAt8A9QHygNrgUaFlvkLMMl/vzOQ6b/fyL98BXyn8P0GKBvGXK2Ayv77dxXk8j/+r4fbqwfwXBHPPR341v9vDf/9GuHKVWj5e4Apod5e/te+ArgE+OQY89sC7wIG/BFYGertVcJczQrWh+9iMysD5m0Eanm0vVoCs0/2PRDsXIWW/TOwKNTbC6gLXOK/XxX4qoj/jyF9f0XaCL0kF9PoAPzTf38mcJWZmX/6DOdcrnNuA7De/3phyeWcW+ycO+B/uALfWSlDrbQXHwl0HTDfObfbOfcjMB9o7VGusF0UxTn3HrD7OIt0AF52PiuA6mZWl9Bur2JzOeeW+9cL4Xt/lWR7HcvJvDeDnSss7y/n3Dbn3Gr//X34zn2VUGixkL6/Iq2gF3UxjcIb5JdlnO/EYT8BNUv43FDmCtQT36dwgYrmOw/8CjO7IUiZSpPrZv+fdzPNrOBUyBGxvazoi6KEanuVxLGyh3J7lVbh95cD/m1mH5pZbw/yXGZma83sXTM73z8tIraXmVXGVxhfD5gc8u1lvlbwxcDKQrNC+v4q0cm5pOTMrCuQDLQImJzonNtiZvWBRWb2sXPumzBFmgVMd87lmlkffH/dXBmmdZdEURdF8XJ7RTQza4WvoF8eMPly//aqA8w3sy/8I9hwWI3v9/VfM2sLvAU0CNO6S+LPwDLnXOBoPqTby8xOxfcB0t85tzdYr1sSkTZCL8nFNH5ZxsxOAU4DdpXwuaHMhZldDQwG2jvncgumO+e2+P/9FliC75M7LLmcc7sCsrwANC3pc0OZK8BRF0UJ4fYqiWNlD+X2KhEzuxDf77CDc25XwfSA7bUDeJPgtRqL5Zzb65z7r//+HKCcmdUiAraX3/HeX0HfXmZWDl8xz3DOvVHEIqF9fwV7x8BJ7lQ4Bd/OgHP4/x0p5xda5m5+vVM0y3//fH69U/RbgrdTtCS5Lsa3E6hBoek1gAr++7WArwnSzqES5qobcP9GYIX7/50wG/z5avjvnx6uXP7lGuLbQWXh2F4B60ji2Dv5rufXO60+CPX2KmGuevj2CzUrNL0KUDXg/nKgdRhz/abg94evMG72b7sSvQdClcs//zR8ffYq4dhe/p/7ZWDMcZYJ6fsraBs3iL+ktvj2Dn8DDPZPexLfqBegIvCa/839AVA/4LmD/c/7EmgT5lwLgO+Bj/y3d/zTmwEf+9/QHwM9w5zraeBT//oXAw0DnnuHfzuuB1LDmcv/+AlgWKHnhXp7TQe2AYfx9Sl7An2Bvv75Boz35/4YSA7T9iou1wvAjwHvr2z/9Pr+bbXW/3seHOZc/QLeXysI+MAp6j0Qrlz+ZXrgO1Ai8Hkh21742mAOWBfwe2obzveXvvovIhIjIq2HLiIiJ0gFXUQkRqigi4jECBV0EZEYoYIuIhIjVNAlbphZdTP7i//+mWY20+tMIsGkwxYlbvjPrzHbOXeBx1FEQkLncpF4Mgw418w+wvcN1POccxeYWQ/gBnzfHGwAjMD37cZuQC7Q1jm328zOxfelkNrAAaCXc+6L8P8YIkVTy0XiyUPAN865JsCgQvMuAG4Cfg+kAQeccxcD/wFu9y+TDtzjnGsKDAQmhCW1SAlphC7is9j5zmG9z8x+wneWSvB9PftC/xn0mgGv+U6/D/jOGyQSMVTQRXxyA+7nBzzOx/f/pAywxz+6F4lIarlIPNmH79JgpeZ857XeYGYd4ZdrQ15UzNNEwkoFXeKG851DfJn/wsL/OIGXSAF6mlnBmfqCdkk1kWDQYYsiIjFCI3QRkRihgi4iEiNU0EVEYoQKuohIjFBBFxGJESroIiIxQgVdRCRG/B+yU8P/8f03ZgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "code", "metadata": { "id": "SinyneEJPMxp" }, "source": [ "# Initial Condition\n", "w=np.zeros(N+1)\n", "#np.zeros(N+1)\n", "w[0]=IC" ], "execution_count": 6, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "OHbrgQYFPMxr" }, "source": [ "\n", "\n", "## 2-step Adams Bashforth (Predictor)\n", "\n", "The 2-step Adams Bashforth difference equation is\n", "\\begin{equation}w^{0}_{i+1} = w_{i} + \\frac{h}{2}(3f(t_i,w_i)-f(t_{i-1},w_{i-1})) \\end{equation}\n", "\\begin{equation}w^{0}_{i+1} = w_{i} + \\frac{h}{2}(3(t_i-w_i)-(t_{i-1}-w_{i-1}))\\end{equation}\n", "## 1-step Adams Moulton (Corrector)\n", "\\begin{equation}w^{1}_{i+1} = w_{i} + \\frac{h}{2}(f(t_{i+1},w^{0}_{i+1})+f(t_{i},w_{i})) \\end{equation}\n", "\n", "\\begin{equation} w^{1}_{i+1} = w_{i} + \\frac{h}{2}((t_{i+1}-w^0_{i+1})+(t_{i}-w_{i})) \\end{equation}\n", "\n", "\n", "For $i=0$ the system of difference equation is:\n", "\\begin{equation}w^{0}_{1} = w_{0} + \\frac{h}{2}(3(t_0-w_0)-(t_{-1}-w_{-1})) \\end{equation}\n", "\\begin{equation}w_{1} = w_{0} + \\frac{h}{2}((t_{1}-w^{0}_{1})+(t_{0}-w_{0})). \\end{equation}\n", "this is not solvable as $w_{-1}$ is unknown.\n", "\n", "For $i=1$ the difference equation is:\n", "\\begin{equation}w^{0}_{2} = w_{1} + \\frac{h}{2}(3(t_1-w_1)-(t_{0}-w_{0})) \\end{equation}\n", "\\begin{equation}w_{2} = w_{1} + \\frac{h}{2}((t_{2}-w^{0}_{2})+(t_{1}-w_{1})). \\end{equation}\n", "this is not solvable as $w_{1}$ is unknown. $w_1$ can be approximated using a one step method. Here, as the exact solution is known,\n", "\\begin{equation}w_1=2e^{-t_1}+t_1-1.$$\n" ] }, { "cell_type": "code", "metadata": { "id": "oMrAer0BPMxr" }, "source": [ "### Initial conditions\n", "w=np.zeros(len(t))\n", "w0=np.zeros(len(t))\n", "w[0]=IC\n", "w[1]=y[1]\n", "w0[0]=IC\n", "w0[1]=y[1]" ], "execution_count": 7, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "Bw_tlt_VPMxs" }, "source": [ "### Loop" ] }, { "cell_type": "code", "metadata": { "id": "FNQoj6TlPMxt" }, "source": [ "for k in range (1,N):\n", " w0[k+1]=(w[k]+h/2.0*(3*myfun_ty(t[k],w[k])-myfun_ty(t[k-1],w[k-1]))) \n", " w[k+1]=(w[k]+h/2.0*(myfun_ty(t[k+1],w0[k+1])+myfun_ty(t[k],w[k]))) " ], "execution_count": 8, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "5rEWmDdZPMxu" }, "source": [ "### Plotting solution" ] }, { "cell_type": "code", "metadata": { "id": "FfQsR7SIPMxv" }, "source": [ "def plotting(t,w,w0,y):\n", " fig = plt.figure(figsize=(10,4))\n", " plt.plot(t,y, 'o-',color='black',label='Exact')\n", " plt.plot(t,w0,'v:',color='blue',label='Adams-Bashforth Predictor')\n", " plt.plot(t,w,'^:',color='red',label='Adams-Moulton Corrector')\n", " plt.xlabel('time')\n", " plt.legend()\n", " plt.show " ], "execution_count": 9, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "NLHumr0SPMxw" }, "source": [ "The plot below shows the exact solution (black) and the Adams Predictor Corrector approximation (red) of the intial value problem" ] }, { "cell_type": "code", "metadata": { "id": "3nDuj6A0PMxw", "outputId": "342d4470-f82a-444f-efab-df73f76b3a1e", "colab": { "base_uri": "https://localhost:8080/", "height": 0 } }, "source": [ "plotting(t,w,w0,y)" ], "execution_count": 10, "outputs": [ { "output_type": "display_data", "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAEHCAYAAACKtsOUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd1hUR/fA8e+F2MUSe6KCGgsqRcWCDUts2JMYTWzYjSYGTWJJYokt+rOX1/haMWqwxdhi76j4xoYdu8aCXcECSpnfHyMoSllhYUHO53l4dmfv3blnkYTDzNwzhlIKIYQQQgiRMFaWDkAIIYQQIjWTZEoIIYQQIhEkmRJCCCGESARJpoQQQgghEkGSKSGEEEKIRJBkSgghhBAiEd6L7wTDMOYBTYDbSqmyMRxvDowAIoAwwFMptSe+fnPnzq3s7OzeOmAhhBBCiOR26NChu0qpPDEdM+KrM2UYRk3gMfB7LMlUVuCJUkoZhuEILFNKlYovKBcXF3Xw4EGTPoAQQgghhCUZhnFIKeUS07F4p/mUUruB+3Ecf6xeZmRZAKkCKoQQQog0wyxrpgzDaGkYhj/wN9A5jvO6G4Zx0DCMg3fu3DHHpYUQQgghLMosyZRS6q8XU3st0OunYjtvllLKRSnlkidPjNOOQgghhBCpSrwL0N+GUmq3YRhFDcPIrZS6+7bvDw0N5dq1a4SEhJgzLCHMJmPGjBQsWJB06dJZOhQhhBApRKKTKcMwPgIuvFiAXh7IANxLSF/Xrl3DxsYGOzs7DMNIbGhCmJVSinv37nHt2jWKFCli6XCEEEKkEPFO8xmG4Q34AiUNw7hmGEYXwzB6GobR88UpnwInDMPwA/4DtFbx3SIYi5CQEHLlyiWJlEiRDMMgV65cMnIqhBApQLlyYBhvfpUrl/yxxDsypZT6Ip7jY4Gx5gpIEimRksnPpxBCpAyurnDqFDx//vK19OmhatXkj0UqoAshhBAi1Rk8GKxey2KsrfXryU2SqddYW1vj7Owc9TVmzBiz9e3n58f69evN1p8QQgiRVhUoAJ066dEo0I+dOkH+/MkfS6pOphYvXoydnR1WVlbY2dmxePHiRPeZKVMm/Pz8or4GDhxohkg1SaaEEEKIxLt+HWbPhtKl/yQ0VK9jDQ0NpnTpPy0ST6pNphYvXkz37t25cuUKSimuXLlC9+7dzZJQvS4wMJCSJUty5swZAL744gtmz54NwFdffYWLiwtlypRh6NChUe85cOAAVatWxcnJiUqVKhEYGMiQIUNYunQpzs7OLF261OxxCiGEEGnBpEnQp08o/fv/gFLzgHCUmkf//h2SJA+IT7x78yWVmPbmO336NPb29gB4enri5+cX6/v379/Ps2fP3ng9Q4YMVKlSJcb3ODs7M3ny5Djjsra2xsHBIao9aNAgWrduzZYtWxgyZAjffvstXl5ebNy4EYD79+/z/vvvEx4eTt26dZk6dSqlSpWiVKlSLF26lIoVKxIUFETmzJlZtGgRBw8eZPr06XHGIFK2V39OhRBCJJ/gYMiUCZ49g6JF63LjxnYgP7AEaA3cwtbWlsuXL5v92nHtzWfWop3JKaZEKq7XTRU5zfe6evXqsXz5cnr37s3Ro0ejXl+2bBmzZs0iLCyMgIAATp06hWEYFChQgIoVKwKQLVu2RMUkhBBCpHVDhsD69eDjoxMqnUgB3ARqRZ3377//JntsKTaZim8Eyc7OjitXrrzxuq2tLTt37jR7PBEREZw+fZrMmTPz4MEDChYsyKVLlxg/fjwHDhwgZ86ceHh4SA0iIYQQIglUrAgPHoCVVQQDB/4Y63mFCxdOxqi0VLtmatSoUWTOnDnaa5kzZ2bUqFFJcr1JkyZhb2/PH3/8QadOnQgNDSUoKIgsWbKQPXt2bt26xYYNGwAoWbIkAQEBHDhwAIBHjx4RFhaGjY0Njx49SpL4hBBCiHeNvz+sW6efN20KEyY8p0uXDowdO5a6desmax4Ql1SbTLVt25ZZs2Zha2uLYRjY2toya9Ys2rZtm6h+g4ODo5VGGDhwIGfOnGHOnDlMmDCBGjVqULNmTUaOHImTkxPlypWjVKlSfPnll1SrVg2A9OnTs3TpUr755hucnJyoV68eISEh1K5dm1OnTskCdCGEEMIE/frBN9/owpyBgYG4u7uzePFiRo8ezZYtW5IkD0iIFLsAXYiUSn5OhRAi6SgFoaG6btStW3qxubX1ddzd3Tl16hTz5s2jffv2yR7XO7kAXQghhBDvFqWgY0cIC4PFiyFfPjh58iSNGjXi4cOHrF+/nnr16lk6zDdIMiWEEEKIFMEwoHRpCA/X7V27dtGiRQsyZcrE7t27cXZ2tmyAsZBkSgghhBAWtX8/ZMwIzs4QufHI0qVL6dChA8WKFWPDhg3Y2tpaNsg4pNoF6EIIIYRI/cLCoH17vdg80sSJE2nTpg2VK1dmz549cSdSAQHg5gY3byZ9sLGQZEoIIYQQyS48XK+Reu89WL0ali/XNR379u3Ld999R6tWrdi8eTPvv/9+3B2NGAF79uhHC5FkSgghhBDJ6ulTcHeHsWN1u3RpyJIlhNatWzN58mQ8PT1ZsmQJGTNmjLsjf3+YPx8iIvSjhUanJJmKwapVqzAMA39//1jPqVWrFq+XdkgqHh4eFClSBGdnZ0qVKsUvv/ySoH6GDRvG+PHjYzw2depU7O3t36o+h5+fH+vXrzep/1fZ2dnh4OCAo6Mj9evX52YifvhfveaQIUPYunWryfEKIYSwjEyZoEAByJNHt+/fv0/9+vVZsWIFEydOZNKkSVhZxZOibN4MTk4vV6uHh1tsdCrVJlPlyulV/69/lSuX+L69vb2pXr063t7eie/MTMaNG4efnx9+fn4sWLCAS5cumbX/GTNmsGXLFpN32w4LC0tUcrJjxw6OHTuGi4sLo0ePjnZMKUVERMRb9zl8+HA+/vjjWI8nJN6wsLC3jkMIIUTMNm7UtaMMA7y8oEsXuHLlCtWrV+d///sfS5cupW/fvqZ19sEHOoEKDdXt588tNjqVapMpV1dd0OtV6dND1aqJ6/fx48fs2bOHuXPnsmTJkqjXg4ODadOmDfb29rRs2ZLg4OCoY1999RUuLi6UKVOGoUOHRr1uZ2fHoEGDcHZ2xsXFhcOHD9OgQQOKFSvGzJkzAQgICKBmzZo4OztTtmxZfHx84owvcu+/LFmyADqBqFixImXLlqV79+5EFmGdOnUqpUuXxtHRkTZt2kS9/9SpU9SqVYuiRYsydepUAHr27MnFixdp1KgRkyZN4v79+7Ro0QJHR0eqVKnCsWPHAD0K1L59e6pVq0b79u0ZMmQIS5cujVbRPab+41KzZk3Onz/P5cuXKVmyJB06dKBs2bJcvXqVcePGUbFiRRwdHaN9X0eNGkWJEiWoXr06Z86ciXrdw8ODFStWAHDgwAGqVq2Kk5MTlSpVIjAw8I14Tf2cQgghEu/uXfjsM3jlf+f4+flRpUoVAgIC2Lx5M59//nncnUREwB9/6McZM8DaOvpxS41OKaUs8lWhQgX1ulOnTkVru7kpNX++fv78uW4vXKjb588rZWWllF6+pr+srJSaM0cfv3NHn79mjW4HBLxxuRgtWrRIde7cWSmllKurqzp48KBSSqkJEyaoTp06KaWUOnr0qLK2tlYHDhxQSil17949pZRSYWFhys3NTR09elQppZStra2aMWOGUkopT09P5eDgoIKCgtTt27dV3rx5lVJKjR8/Xo0cOTLq/UFBQW/E1LFjR2VnZ6ecnJxUlixZ1KBBg6KORV5bKaXatWun1rz4wAUKFFAhISFKKaUePHiglFJq6NChytXVVYWEhKg7d+6o999/Xz1//jwq1jt37iillPr666/VsGHDlFJKbdu2TTk5OUW9v3z58urp06dKKaXmz5+vevfuHXX9uPp/1avX6t27t+rfv7+6dOmSMgxD+fr6KqWU2rRpk+rWrZuKiIhQ4eHhqnHjxmrXrl3q4MGDqmzZsurJkycqMDBQFStWTI0bNy7q+7R8+XL17NkzVaRIEfXPP/8opZQKDAxUoaGhb8Rr6ud83es/p0IIIWIXHv7yuY+PUk+e6OebN29WWbNmVYUKFVInTpwwrbPVq/Uv/JUrlXJ2jp4ERH45O5v/QyilgIMqlpwm1Y5MFSigK6O+96JSVvr0kD8/5MyZuH69vb2jRnLatGkTNdW3e/du2rVrB4CjoyOOjo5R71m2bBnly5enXLlynDx5klOnTkUda9asGQAODg5UrlwZGxsb8uTJQ4YMGXj48CEVK1Zk/vz5DBs2jOPHj2NjYxNjXJHTfDdv3mTbtm3s27cP0NNllStXxsHBge3bt3Py5MmoGNu2bcuiRYt4772X5cQaN25MhgwZyJ07N3nz5uXWrVtvXGvPnj1RIzJ16tTh3r17BAUFRX2eTJkyxfr9M6V/gNq1a+Ps7ExQUBCDBg0CwNbWlipVqgCwefNmNm/eTLly5Shfvjz+/v6cO3cOHx8fWrZsSebMmcmWLVvU9/dVZ86coUCBAlSsWBGAbNmyRfsemONzCiGEiN+tW3omKXKz4urVIXNm+P3333F3d6do0aLs37+fMmXKmNZh06awaRO0aAFHjsSUSunXk1mKLtq5c+fL5+nSRW9nzgyHDkHRorpGhbW1bufPr4/nzh39/MjX43L//n22b9/O8ePHMQyD8PBwDMNg3Lhxsb7n0qVLjB8/ngMHDpAzZ048PDyipuIAMmTIAICVlVXU88h2WFgYNWvWZPfu3fz99994eHjQr18/bGxsohaZz5kzJ9r1smbNSq1atdizZw/ly5enV69eHDx4kEKFCjFs2LCoa//999/s3r2btWvXMmrUKI4fPx4tHgBra+u3XhMUOb0YG1P737FjB7lz545qP3z4MFrfSikGDRpEjx49or1v8uTJbxVvQsX3OYUQQsTPxkb/vo6klOLXX3/lp59+om7duqxcuZJs2bLF3UlAAHTrpqf1CheG+vWTNugESLUjU6BHpzp1Aisr/WhKwhSXFStW0L59e65cucLly5e5evUqRYoUwcfHh5o1a/LHH38AcOLEiaj1NUFBQWTJkoXs2bNz69YtNmzY8FbXvHLlCvny5aNbt2507dqVw4cP07Jly6jF5i4u0fdUDAsL43//+x/FihWLSpxy587N48ePo9YLRUREcPXqVWrXrs3YsWMJDAzk8ePHJsdUo0aNqIXoO3fuJHfu3DH+sNvY2PDo0aO3+rymatCgAfPmzYuK+/r169y+fZuaNWuyatUqgoODefToEWvXrn3jvSVLliQgIIADBw4A8OjRI8LCwt6I19TPKYQQ4u2sXavXg2fODNu3Q5Mm+vdXr169+Omnn2jXrh3r16837f+5d++Cnx9cvJj0gSdQih6ZMsXgwXDypH5MLG9vbwYMGBDttU8//RRvb28mTpxIp06dsLe3x97engoVKgDg5OREuXLlKFWqFIUKFaJatWpvdc2dO3cybtw40qVLR9asWfn9999jPO+HH35g5MiRPH/+nLp16/LJJ59gGAbdunWjbNmy5M+fP2paKzw8nHbt2hEYGIhSij59+pAjRw6TYxo2bBidO3fG0dGRzJkzs2DBghjPq127NmPGjMHZ2Tlqqs5c6tevz+nTp3F1dQX0iNyiRYsoX748rVu3xsnJibx580Z95lelT5+epUuX8s033xAcHEymTJnYunXrG/Ga+jmFEEKYzs8PmjWDSZPA01PfuffkyRO++OIL1q5dy8CBAxk9ejSGYcTd0c2bepTEwQHOn9f7zaRQhnpx91dyc3FxUa/XaTp9+jT29vYWiUcIU8nPqRBCvEkpnTgB/P03NGig1zXfuXOHJk2acPDgQaZNm0avXr3i7+yff6B2bViwQN8CmAIYhnFIKeUS07FUPc0nhBBCCMs7cwYqVoTTp3W7cWOdSF24cIGqVaty7NgxVq5caVoiBboYZ7duULNm0gVtRvEmU4ZhzDMM47ZhGCdiOd7WMIxjhmEcNwxjn2EYTuYPUwghhBApVebM+mawFzdEA/DPP//g6urKgwcP2L59O82bN4+7E6Vg1iy910yGDDB5MuTNm7SBm4kpI1NeQMM4jl8C3JRSDsAIYJYZ4hJCCCFECqYUrF+vHwsVgsOHoXJlfWzdunXUrl2brFmzsm/fvqj1r3E6fBh69tSl0VOZeJMppdRu4H4cx/cppR68aO4HCpopNiGEEEKkUMuW6em8jRt1O3IrvVmzZtG8eXNKly6Nr68vJUqUiLujyLXbFSqAry989VXSBZ1EzL1mqgsQa20AwzC6G4Zx0DCMg3fu3DHzpYUQQgiR1CJzn1atdELVsGHk64rBgwfTo0cPGjZsyI4dO8iXL1/cnd25A3Xr6gXnoIe24rvLLwUyWzJlGEZtdDI1ILZzlFKzlFIuSimXPJFbRQshhBAiVfD11WvC79/XI1GtWuncJzQ0lE6dOjFy5Ei6dOnC6tWryZo1a/wdhofD7ds6qUrFzJJMGYbhCMwBmiul7pmjT0tatWoVhmHg7+8f6zm1atXi9dIOSWXnzp0YhhGtGrqfnx+GYTB+/PgE99mkSZOo55Hb05jDzZs3adOmDcWKFaNChQq4u7tz9uxZs/UfHy8vL27cuJFs1xNCiLTkwQOdTEV69OgRTZo0YcGCBQwbNozZs2fHuIVXNBcv6iGu/Pnh6FE9X5iKJTqZMgyjMLASaK+USr7fmJECAsDNTRf3MhNvb2+qV68etS9fSlC2bFmWLVsW1fb29sbJyTw3TpozmVJK0bJlS2rVqsWFCxc4dOgQv/76a6x79L3u9e1n3na7G0hYMhUeHv7W1xFCiLQgLAx279bPXV117vPRR7odEBCAm5sb27ZtY+7cuQwdOjT+YpznzoGjI0QOBlhbJ13wycSU0gjegC9Q0jCMa4ZhdDEMo6dhGD1fnDIEyAXMMAzDzzCM5BmuiTRiBOzZox/N4PHjx+zZs4e5c+eyZMmSqNeDg4Np06YN9vb2tGzZkuDg4KhjX331FS4uLpQpU4ahQ4dGvW5nZ8egQYNwdnbGxcWFw4cP06BBA4oVK8bMmTMB/YNYs2ZNnJ2dKVu2LD4+PjHGZWtrS0hICLdu3UIpxcaNG2nUqFHUcT8/P6pUqYKjoyMtW7bkwQN9T8CrI2h3797Fzs4uWr+XL19m5syZTJo0CWdnZ3x8fLh8+TJ16tTB0dGRunXr8u+//wLg4eFBnz59qFq1KkWLFo3avuZVO3bsIF26dPTs2TPqNScnJ2rUqIFSih9++IGyZcvi4ODA0qVLAZ3M1ahRg2bNmlG6dOk32uHh4fzwww9UrFgRR0dH/vvf/0b1PXbsWBwcHHBycmLgwIGsWLGCgwcP0rZtW5ydnQkODmbbtm2UK1cOBwcHOnfuzLNnz6L+fQYMGED58uVZvnx5bD8SQgiRpo0eDXXq6BwIXuY+/v7+uLq6cvbsWdatW0fnzp1N6/Cjj+DHH6Fdu6QJ2BKUUhb5qlChgnrdqVOnor/g5qbU/Pn6+fPnur1woW4/eaJUlSpKpUun94nOmFEpV1el/vxTH79zR5+/Zo1uBwS8cb2YLFq0SHXu3FkppZSrq6s6ePCgUkqpCRMmqE6dOimllDp69KiytrZWBw4cUEopde/ePaWUUmFhYcrNzU0dPXpUKaWUra2tmjFjhlJKKU9PT+Xg4KCCgoLU7du3Vd68eZVSSo0fP16NHDky6v1BQUFvxLRjxw7VuHFjNWXKFDVt2jS1Z88e5eHhoYYOHarGjRunlFLKwcFB7dy5Uyml1ODBg9W333774lvoFhXnnTt3lK2tbbQ+lVLR+lFKqSZNmigvLy+llFJz585VzZs3V0op1bFjR/XZZ5+p8PBwdfLkSVWsWLE3Yp0yZYry9PSM8Xu7YsUK9fHHH6uwsDB18+ZNVahQIXXjxg21Y8cOlTlzZnXx4sWo2F5t//e//1UjRoxQSikVEhKiKlSooC5evKjWr1+vXF1d1ZMnT6L9O7z6mYODg1XBggXVmTNnlFJKtW/fXk2aNCnq32fs2LExxhqXN35OhRDiHRYYqNTSpdFf8/HxUTlz5lT58uWL+j0Zp4gIpaZNU+rGjaQJMhkAB1UsOU3qroB+5crL2woiInQ7kby9vWnTpg0Abdq0iZrq2717N+1eZNGOjo44OjpGvWfZsmWUL1+ecuXKcfLkSU6dOhV1rFmzZgA4ODhQuXJlbGxsyJMnDxkyZODhw4dUrFiR+fPnM2zYMI4fP46NjU2ssX3++ecsX74cb29vvvjii6jXAwMDefjwIW5ubgB07NiR3ZFjsgng6+vLl19+CUD79u3Zs2dP1LEWLVpgZWVF6dKlTZ66i7Rnzx6++OILrK2tyZcvH25ublGbEVeqVIkiRYpEnftqe/Pmzfz+++84OztTuXJl7t27x7lz59i6dSudOnUi84styd9///03rnnmzBmKFCkSdWvu69+b1q1bv9VnEEKItGDtWvj0Uz3Fly0bfP75y2N//vknH3/8MXny5MHX1zdqr9o4/fsvDBwIL2Zl3jUpe6PjnTtfPk+XLno7MFCvgotcU/P8uW5XrarbuXNHPz9//ngvd//+fbZv387x48cxDIPw8HAMw2DcuHGxvufSpUuMHz+eAwcOkDNnTjw8PAgJCYk6niFDBgCsrKyinke2w8LCqFmzJrt37+bvv//Gw8ODfv36YWNjwy+//AIQbdF5/vz5SZcuHVu2bGHKlCkmrXN67733iIiIAIgWV0K9+hlUDPs6lilTJsbpv/hkyZIl1rZSimnTptGgQYNo52zatOmtrxPfdYUQQsDdu3p84uFD/es00tSpU/H09KRKlSqsWbOG3K8ejElEhL7tz9YWDhyAkiWTNnALSb0jUyNG6H+kV4WHJ2rt1IoVK2jfvj1Xrlzh8uXLXL16lSJFiuDj40PNmjX5448/ADhx4gTHjh0DICgoiCxZspA9e3Zu3brFhg2xltmK0ZUrV8iXLx/dunWja9euHD58mJYtW+Ln54efnx8uLtH3VBw+fDhjx47F+pUFe9mzZydnzpxR660WLlwYNUplZ2fHoUOHoj5fTGxsbHj06FFUu2rVqlHrxRYvXkyNGjVM/jx16tTh2bNnzJr1shD+sWPH8PHxoUaNGixdupTw8HDu3LnD7t27qVSpUrx9NmjQgN9++43Q0FAAzp49y5MnT6hXrx7z58/n6dOngE6GX/88JUuW5PLly5w/f/6N740QQoiXnj4FPz/9vFMnXQYhMleKiIjghx9+4Ntvv6V58+Zs27Yt/kTq4UO9WXHk7x57+5eVPd8xKXtkKi6+vno06lXPn0Mi7krz9vZmwIDoZbI+/fRTvL29mThxIp06dcLe3h57e/uoYU0nJyfKlStHqVKlKFSoENWqVXura+7cuZNx48aRLl06smbNyu+//x7n+VUjR95es2DBAnr27MnTp08pWrQo8+fPB+D777/n888/Z9asWTSO5dbTpk2b8tlnn7F69WqmTZvGtGnT6NSpE+PGjSNPnjxRfZnCMAz++usvPD09GTt2LBkzZsTOzo7JkydTvXp1fH19cXJywjAM/u///o/8+fPHWYICoGvXrly+fJny5cujlCJPnjysWrWKhg0bRiWc6dOnx93dndGjR+Ph4UHPnj3JlCkTvr6+zJ8/n1atWhEWFkbFihWjLY4XQgihde8OmzbpqgU2NnpCCODZs2d4eHiwZMkSevfuzZQpU6L9QR+rDBn0bscxzGK8a4yYpmqSg4uLi3q9TtPp06ext7e3SDxCmEp+ToUQ76ILF8DfP3rJp4cPH9KyZUt27tzJ2LFj+eGHH+IvfeDvD0WK6GRKqVRZ0TwmhmEcUkq5xHTs3RxvE0IIIUS85s2D/v3182LFoidSV69epUaNGuzdu5fFixfTv3//+BOpu3ehShX4/nvdfkcSqfik3mk+IYQQQiTK8eNw8qReJZM+/auvH6dRo0Y8evSIjRs3UqdOHdM6zJ0bJk/W++2lISkumVJKxZ/5CmEhlpoWF0IIc7l3Dx49Ajs7iLxZ/dXdX7Zv307Lli2xsbHBx8cnWimgWE2bBjVqgLMzeHgkRdgpWoqa5suYMSP37t2TX1giRVJKce/ePTJmzGjpUIQQIkGUAnd3XUMqIkInUa8mUn/88QcNGzakUKFC+Pr6mpZIBQbC//0fvHIXd1qTokamChYsyLVr17iTynePFu+ujBkzUrBgQUuHIYQQCWIYMHYsZMwYvUqBUopx48YxYMAA3NzcWLVqFTly5Ii7s9BQnYllz67vpP/gg6QNPgVLUclUunTpolXBFkIIIUTiKKX31ytQADp3hlq1oh8PDw/H09OT6dOn06ZNG7y8vKIVaI7R06fQtKleG/Xjj1CoUJLFnxqkqGk+IYQQQphXeDjs2gV79755LDg4mFatWjF9+nS+//57Fi9eHH8iBXpoq3DhNJ9ERUpRI1NCCCGEMI+rV/UMXLZssGoVZMoU/fi9e/do2rQp+/fvZ8qUKfTp0yf+Tv39IVcuyJMH3qKg87tORqaEEEKId8yTJ7rc01df6XbmzNFLPl26dImqVaty+PBhli9fbloiFRICH38MXbsmTdCpmIxMCSGEEO+YLFn0QnOXGOp1Hzp0iMaNG/P8+XO2bt1K9erVTes0Y0aYMwdKlDBvsO8AGZkSQggh3gGhodC378u1Ue3aQalS0c/ZuHEjbm5uZMyYkb1795qWSM2cCWvX6ucNG0LRouYN/B0gyZQQQgjxDnjyBNatg+3bYz4+f/58mjRpQokSJfD19TVtj9HQUJg7FxYtMm+w7xiZ5hNCCCFSsYsXdTXzHDng0CG94PxVSilGjBjB0KFDqV+/PitWrMDGxibuTkND9WO6dLB5M8R3fhonI1NCCCFEKnXpEjg5wZgxuv16IhUWFkb37t0ZOnQoHTt2ZN26dfEnUuHh0Lw5dOumi1TlzBm9TLp4g3x3hBBCiFTKzg4GD4a2bd889vjxY1q3bs369ev5+eefGT58uGl731pbQ/XqkDdv9FsARawMS+2D5+Liog4ePGiRawshhBCp1ZMn8N138NNPsdfMvHXrFk2aNOHw4cPMmDGDHmVLfWgAACAASURBVD16xN/x+fN6es+UtVRpkGEYh5RSMdwfKdN8QgghRKpy7RosXaqrmsfk7NmzVK1alZMnT7Jq1SrTEqmICL37cbt2empPvBWZ5hNCCCFSgX//1Tu4lCwJFy7A+++/ec7+/ftp0qQJhmGwc+dOKlWqZFrnVlawYMGb1T2FSWRkSgghhEjhdu+Gjz7S28JAzInU6tWrqVOnDjly5MDX19e0RMrLC6ZM0c+dnaUgZwJJMiWEEEKkcFWq6HVStWrFfPy3337jk08+wcHBAV9fXz766KP4O1UK/v4bNmzQ03wiweJNpgzDmGcYxm3DME7EcryUYRi+hmE8Mwzje/OHKIQQQqQ9d+/C11/D06eQPj38+quuJfUqpRSDBg2iV69euLu7s337dvLkyRN3x+Hh8OiRns5buBDWrNHTfCLBTPnueQEN4zh+H+gDjDdHQEIIIYTQBTjnzYMDB2I+/vz5czp06MCYMWPo0aMHf/31F1myZIm/47ZtoXFjCAvT++2lT2/ewNOgeBegK6V2G4ZhF8fx28BtwzAamzEuIYQQIk0KCIACBaBBA7h8WZd7el1QUBCffvopW7duZeTIkfz444+m1ZACaNkSbt+WQpxmJN9JIYQQIoVYuBB69oT9+8HBIeZE6vr167i7u3Pq1Cm8vLzo2LFj/B1fuaIzMzc3aN3a7HGndcmaTBmG0R3oDlC4cOHkvLQQQgiR4jVoAF99BcWLx3z85MmTNGrUiAcPHvD3339Tv3590zru0QNOn4Zz52RaLwkk64ozpdQspZSLUsol3gVyQgghRBpw5Qr8/LO+uS5vXhg/Xi9let2uXbuoXr06oaGh+Pj4mJ5IAcydC+vXSyKVRGT5vhBCCGFBf/0F06frQpyxWbZsGfXr1yd//vzs378fZ2fn+Dv29oY+fXSW9uGHUKaM+YIW0ZhSGsEb8AVKGoZxzTCMLoZh9DQMo+eL4/kNw7gG9AN+fnFOtrj6FEIIIdK6u3f147ffwsmTuihnTCZNmkTr1q2pVKkSe/fuxdbW1rQLHD8Ox45BSIh5Ahaxko2OhRBCiGQ2bJieeTt8GGJb9RIREcF3333H5MmT+fTTT1m0aBEZY5r/i/4muHMH8uXTI1LPn0OGDGaPPy2Ka6NjuZtPCCGESGYtWug8J6ZtYQBCQkLo0KEDy5cvp0+fPkycOBFra+v4O+7dG7ZsgSNHwMZGEqlkIsmUEEIIkQyOH4e9e3XpA2dn/RWT+/fv06JFC3x8fBg/fjz9+vUzvYZUx476VsCsWc0XuIiXJFNCCCFEMpgyBTZt0gXIbWxiPufKlSs0atSICxcu4O3tTZs2beLv+Pp1naV9/rnexK9KFfMGLuIlyZQQQgiRRCIi4PFjyJYNpk2Dhw9jT6T8/Pxwd3fn6dOnbNq0iVqx7Wr8umHDYMUK+Pjj2OcNRZKS0ghCCCFEEunUCRo1gtBQyJRJbxMTky1btlCzZk2sra3Zu3evaYlU5A1kkyeDj48kUhYkyZQQQgiRRJo0gU8/jXsbvIULF+Lu7o6dnR379++njCn1oFauhObN9Sr2LFmgbFnzBS3emkzzCSGEEGa0a5ee2mvcGFq1iv08pRRjxozhxx9/pE6dOqxcuZLs2bObdpGgIF2o6skTqWqeAkidKSGEEMJMlILq1fW03v79YBXL/E9YWBjffPMNM2fO5Msvv2T+/Pmkjy8pUgouXYKiRXU7PBxMKZcgzCKuOlMyzSeEEEIk0vPn8OwZGIZeC75lS+yJ1NOnT/nkk0+YOXMmAwcOZOHChfEnUgBDhkCFCvruPZBEKgWRaT4hhBAiEUJDoX59KFZMVzWPbZE5wJ07d2jatCn//PMP06dPp3fv3qZfqGtXyJkTPvgg8UELs5KRKSGEECIR0qWDevWgbt24z7tw4QJVq1bl6NGjrFy50rRE6tYtfbeeUmBrC/366eEvkaLIyJQQQgiRAH/+CSVL6hvpfvop7nMPHDhA48aNCQ8PZ9u2bVStWtW0i8yeDb/+Cs2avVwrJVIcGZkSQggh3tKTJ/DNNzB6dPzn/v3339SqVYusWbOyb98+0xKpiAj9+OOPejdkSaRSNEmmhBBCCBOFhOgZtyxZYPt2mD8/7vNnz55Ns2bNsLe3Z9++fZQsWTL+i6xfD5Uqwf37ehW7Ke8RFiXJlBBCCGGCe/d0jjN1qm6XKgUZMsR8rlKKIUOG0L17d+rXr8/OnTvJnz+/aRfKkkV3HBpqnsBFkpNkSgghhDBBzpzg4gL29nGfFxoaSufOnRkxYgSdO3dmzZo1ZM2aNe43KQVHj+rnbm6wZw/ky2eewEWSk2RKCCGEiIOXlx6VsrKCefN0GYTYPHr0iKZNm+Ll5cWwYcOYM2cO6dKli/8i06frTC0yoZI79lIVuZtPCCGEiMXFi9CjB1y+DMOGxX3uzZs3ady4MUePHmXOnDl06dLF9At17KgfHR0TGqqwIBmZEkIIIV7z/Ll+LFoU9u6FwYPjPt/f3x9XV1f8/f1Zs2aNaYnU3bv6br2wMMiWTd8eKCNSqZIkU0IIIcQrzp6F0qVh61bddnGJe+eWvXv3Uq1aNZ4+fcquXbtwd3c37UIbN8KkSeDnl/ighUVJMiWEEEK84oMPoEQJveA8PitXrqRu3brkypULX19fXFxi3Ac3usi79Nq1gzNndLYmUjVJpoQQQqR5ERF6cXloKGTNqks9VagQ93umTZvGZ599Rrly5di3bx9FTSmsuXOnrht19qxuFy6c6NiF5UkyJYQQIs3bvh26dIHly+M/NyIigv79+9OnTx+aN2/Otm3byJ07t2kX+vBDvRArS5bEBSxSFEmmhBBCpFlhYfrx449h1y744ou4z3/27Bnt2rVj3Lhx9OrVixUrVpA5c+a436QU+Pjo58WL68VYH36Y+OBFiiHJlBBCiDRp3z5dgPPcOd2uWTPum+kePnxIo0aN8Pb2ZsyYMUyfPh3ruFamR1q+XHe+aZN5AhcpjtSZEkIIkSYVKKAXm5tSjeDq1au4u7tz5swZFi5cSLt27Uy/0CefwJw5UK9ewoMVKVq8I1OGYcwzDOO2YRgnYjluGIYx1TCM84ZhHDMMo7z5w3w75crp/zgKGAHsMtzIb9zEMPTrQggh0q7nz2HpUv28SBE9tffRR3G/5/jx47i6unLlyhU2bNhgWiL14AH07AmBgfDee3pBlpVMBr2rTPmX9QIaxnG8EVD8xVd34LfEh5U4uXKdBUIYzAiqs4fBjABCyJXrjKVDE0IIYUGzZ0ObNnDwoGnn79ixg+rVq6OUwsfHh7p165r2xmPHYNEi+OefhAcrUo14kyml1G7gfhynNAd+V9p+IIdhGAXMFWBC+Pu3JT836MR8rImgE/PJxw3OnHmLYVkhhBDvjIgI/dizp17/bUppJ29vbxo0aEDBggXx9fXFyckp/jc9faof3dzg0iWZ2ksjzDHm+CFw9ZX2tRevvcEwjO6GYRw0DOPgnTt3zHDpmN24cYjBeJAevR+AFeEMxoPr1w8l2TWFEEKkTH/9BZUqQVCQrmQe3+CSUopx48bx5ZdfUrVqVfbs2UNhU+pBHTyoyx7s2qXbefIkPniRKiTrBK5SapZSykUp5ZInCX/IKnzwAZ3wwRr9p0hGntMJH8p/8EGSXVMIIUTKlCOHLsQZEhL/ueHh4fTp04f+/fvz+eefs2nTJnKaUgoddCJVrRrY2SUqXpH6mCOZug4UeqVd8MVrFrO4VClevznjPWBJgCLg7A1LhCSEECIZPXr0shJB7dqwYwfkzRv3e4KDg2nVqhXTp0+nX79+eHt7kyFDhvgvtn49hIfD++/Dn3+CrW3iP4BIVcyRTK0BOry4q68KEKiUCjBDvwlW4t49Mr72WnqgYMQdmjZqyunTpy0RlhBCiGQyaBC0bAm3bul2fOUP7t27x8cff8yqVauYNGkSEyZMwMqUu+/27IHGjcHLK9Exi9TLlNII3oAvUNIwjGuGYXQxDKOnYRg9X5yyHrgInAdmA72SLFpTHTmiK86+9uWzZh3/PrpKJZdKeHWabOkohRBCJFJkKZzXv3bvhg0bIF+++Pu4dOkS1apV49ChQyxduhRPT0/TA6heXY9GeXgk+DOI1M+Uu/m+UEoVUEqlU0oVVErNVUrNVErNfHFcKaV6K6WKKaUclFIm3nCa/Oo1rc/hw4f5OlM9PLz6MrT+d4RF7iUghBAi1XF1hfTpo7+WPj3UqKFvqIvP4cOHcXV15fbt22zZsoVWrVrF/6ZHj6BdO323HuiinKZUQhfvrDRXQaxgwYL8eGYBY0p3Y/iWidSrV49bN29ZOiwhhBAJMHjwm7Uwra316/HZuHEjNWvWJEOGDOzdu5caNWqYdtHr12HzZjgkd4gLLc0lUwA2uWwYeHIWCxYs4Mq+i1z7oA5rBntbOiwhhBBv4cED+PprPQoVOTqVPj106gT588f93vnz59OkSRM++ugjfH19sbe3N+2CAKVKwYUL8NlnifsA4p2RJpOpSB06dGDm2OmkI5yffx3Gb7/9hlLK0mEJIYSIw717+jF7dggIgI8/fjk6Fd+olFKK4cOH07lzZ+rUqcPu3bv5wJSyOf7+UKIELFyo2zY2ifsQ4p2SppMpgPqeTfng5h4K1i9Gr169GFa2K/evx1XwXQghhKX07QsVKkBoqE6g9u6F/v31aJSVVdyjUmFhYfTo0YOhQ4fSoUMH1q1bR7Zs2Uy7cNGi+vZAV1fzfRjxzkjzyRRA7ry5WbduHb98+RNDTnkxu1RHLkUuLBRCCGExSukyTk+e6HaTJtCnz8vtYSJLHgwerG+si21U6smTJ7Ro0YLZs2fz448/4uXlRfrXV67H5M8/4fFjPX84a1b8uyKLNEmSqResrKwYsngkS7pOZYLV/6hQoQLr/lpn6bCEECJNO3JEl3FasEC369aFfv3g9VqaBQroXVxiGpW6ffs2tWrVYsOGDfz222+MGjUKI77CUwDnzkHr1jBpUuI/iHinGZZaI+Ti4qIOmrptdzK7cOECrZu3ZsbJYI6WqEaX0zNNK94mhBAi0by9ITBQb0oMsG4dNGgA6dK9fV/nzp2jYcOGBAQEsHTpUpo2bfp2HezcqbeIScjFxTvFMIxDSqkYt8iWDCEGxYoVY9PGTVzLWpg1Z/1o1qwZDyLv4hBCCGF2r/5dv3IlLF788rUmTRKWy+zfvx9XV1eCgoLYsWOHaYnUkyfw+efg66vbtWpJIiXiJclULHIVzEXLoPU0+o8HmzdvxrNIGzb932pLhyWEEO+cHTt0tYGAFxuRzZ6tp+xMmYmLzZo1a6hTpw45cuRg3759VK5c2bQ3Pn0Kx46BbDsm3oIkU3EwDINevXqxbdM2fgo8izFgBIsWLbJ0WEIIkeo9fgy3b+vnhQvrtU73X9xInSPHm4U438Zvv/1Gy5YtKVu2LPv27aN48eLxv+nmTT0UlicPHD0KnTsnPACR5kgyZYIatWsQvmMdUyvkpX379nzd5Wse339s6bCEECJVev4c7O1h4EDdLlZMj0SVKZOw/hYvXoydnR1WVlZkz56dXr160ahRI3bs2EHevHnj7yAgAJydYfhw3X59dbsQ8ZBkykT2tcqwav8avvvuO1znHeFY3ppcuXDF0mEJIUSqcPcuRA7sp08Pw4ZB9+6J73fx4sV0796dK1euoJQiKCiI9957j9atW5MlSxbTOsmfX692//zzxAck0iS5my8BFjQZxcmN21iQ6yTLli3DzZTdNIUQIg375Rc98HP5MhQqZL5+7ezsuHLlzT9sbW1tuXz5ctxvXrpU36lXsKD5AhLvLLmbz8w6rvsJj2PTyZkzJ9/X7sXcsr2ICI+wdFhCCJFi3L8PvXvrCuWgC22eOGHeROrGjRsxJlIA//77b9xvvntXD42NHm2+gESaJclUApUuXZp//vmHr7MVo97J1bRv1o7Hj2UdlRAibQsN1Y8ZM8KaNbroJkDOnHqdlDkEBwczatQoSpQoEes5hQsXjvlA5GxM7tx6oZYU5BRmIMlUImTLlo3291bh/U1vlmxcSkWXivgu87F0WEIIYREDBoCbm85XMmfWBcS//tp8/SulWLFiBfb29vz888/Ur1+fiRMnkjlz5mjnZc6cmVGjRr3ZQUiIXhe1ZIluOzvLYnNhFpJMJZKVtRUDpv7Ili1baHapIGVaN2bhkHmWDksIIZKFnx+Eh+vnZcpAjRr6bj3Qo1PmcuTIEWrVqkWrVq3Inj0727dvZ+XKlfTt25dZs2Zha2uLYRjY2toya9Ys2rZtG3NHd+++rMkghJnIAnQzOrLmANvbDef7R+sYMGAAI0eO5L333rN0WEIIkSR279YjUUuW6C3sksKtW7f4+eefmTt3Lrly5WLkyJF07doVa2tr0zu5fh3efx8yZdKZ39u8V4gXZAF6MinXrCJf31lBjx49mD92LqttGnB2r7+lwxJCCLNQCrZtg7Vrdbt6dZg+HRo2NP+1nj17xrhx4yhevDheXl707duXc+fO0aNHj7dLpB4/BldX+Oor3ZZESiQBGZlKIrNajOLT1RP47P2SjN0whUqVKlk6JCGESBSloGpVXZ088i49819DsWaNrul34cIFmjRpwoQJE+JcbB6vOXOgUiVwdDRfoCLNkZEpC+i+6ifObFzLJZsAatSowX+6jJPyCUKIVGfPHqhXT+//axjg7a1Hp5LCiRMnqFevHi1atCB9+vRs3LiRtWvXJiyRWr4cDh3Sz7t2lURKJClJppJQ1QbVOHToEO2LtaD3vP6Mte9KcHCwpcMSQog4hYXp/X5Bj0JduaKLbQLY2Zl3YTnA3bt36dWrF05OThw+fJhp06Zx9OhRGjRokLAOQ0L0rYW//mreQIWIhSRTSSxXrlzMOLSQGfY9GH7uD6pXrx5/VV4hhLCQJ090PajIPKRqVfD3T/i+eXEJDQ1l8uTJFC9enFmzZtGrVy/OnTvH119/Tbp06d6us4AAqFlTP2bMCNu3w8KF5g9aiBhIMpUM0mdKT69TM1m2Zjn/nvuXk0Vb491puqXDEkIIAIKDYedO/TxLFn1nnqvry+NWSfCbYsOGDTg4ONC3b18qVarE0aNHmTZtGrly5UpYh7/8Aj4+0KyZbtvZ6bv3hEgGkkwlo6ZNm7J+4VoKEcifXosZOXIkERGyjkoIYVmDBuk78u7c0e2RI8HdPWmu5e/vj7u7O+7u7kRERLB27Vo2btxImcQMfQUEwIIF+rmfH9y8aZ5ghTCRJFPJrGLzKtje+R8ZvizK4MGD6e3Unn+Pxby3lBBCJIWgIBgxAs6e1W1PT9i8GfLkSbprPnjwAE9PTxwcHNi3bx8TJkzgxIkTNGnSBMMwEtbp1au6ovn330PkH6ZWVvrDCZGMTEqmDMNoaBjGGcMwzhuGMTCG47aGYWwzDOOYYRg7DcOQLbjjkD1XdhYtWsSk4ZMYeWID+yt04Pjx45YOSwjxjoushBMcrNdEbdyo23Z2erlRUggLC2PGjBkUL16cadOm0aVLF86dO0e/fv1Inz594jq3toYdO2DFipdl158/h/nzZXRKJKt4kynDMKyB/wCNgNLAF4ZhlH7ttPHA70opR2A4ILdQxMMwDDwHe7L7u0n8mv0BlStXZtHviywdlhDiHfXLL9C+vX6eLx9cugR9+iTtNbdu3Uq5cuXo3bs3jo6OHDlyhJkzZ5InMUNgU6dC5876+QcfwCefvHlOeLiMTolkZcrIVCXgvFLqolLqObAEaP7aOaWB7S+e74jhuIhFy/Ed2XBiMxXKV+B5x5kszvMZz0KeWTosIcQ7ICDg5XNra3jvvZf76OXLl3TXPX/+PM2bN6devXo8ffqUlStXsm3bNhwTWuvp1eLSgYF6f73Ikah//nn5PNLz57BvX8KuJUQCmJJMfQhcfaV97cVrrzoKRP550BKwMQzjjVsyDMPobhjGQcMwDt6JXOkoyJ8/P1u3bCV93gJcuHuXuh/X5caNG5YOSwiRim3ZAoUL6xvcAH7+Gby8knY3lcDAQH744QdKly7N9u3bGTNmDCdPnqRly5YJXxd14QK4uLz8ID/9BGvWQOQU4ZEjOtl6/evIEfN8KCFMYK4F6N8DboZhHAHcgOtA+OsnKaVmKaVclFIuiRrmfQdlyJSBdreWU3xxD44cOULr0s3589t5lg5LCJGKHD8OBw7o59WqwXffQbFiSX/d8PBwZs+eTYkSJZgwYQLt27fn3LlzDBgwgIwJrfAZGqofCxTQdaNerSIqRAoT7958hmG4AsOUUg1etAcBKKViXBdlGEZWwF8pFeci9Hd9b77EOHHiBNecumIfcY2Vv/bDc0DfhP9VJ4RIEyIioEQJPRq1fXv855vLrl278PT0xM/Pj2rVqjFlyhQqVKiQuE6HDdOr4/ftk+RJpBiJ3ZvvAFDcMIwihmGkB9oAa167QG7DMCL7GgTIkEoilC1bljJHlzG8ohv9Bn3Hl198yd2rdy0dlhAihTl0CLp31+ugrKxgyRJ9Y1tyuHTpEq1ataJWrVrcv3+fJUuW4OPjk/BE6smTl+UNSpTQGxM/k/WjInWIN5lSSoUBXwObgNPAMqXUScMwhhuG8aLULLWAM4ZhnAXyAaOSKN40o1DZwszev5DRo0eTb+kNbtlW5dCmA5YOSwhhYUq9XER+6RKsWgXnz+u2iwu8/37SXv/Ro0f8+OOP2Nvbs379eoYPH46/vz+tW7dO+Aj6xYtQvDgsXqzbX36p79qTCuYilYh3mi+pyDSf6ZZ0nUHw/GV8m/kwCxctpHlzuVlSiLTowQNo1Ag8PKBnTz2QExICmTMn/bUjIiJYuHAhgwYNIiAggHbt2vHrr79SsGAiygreuqVvK4yI0B+oSxeoXNl8QQthRomd5hMW1mZOL+pcXEDJUiXp0qITs4p25nnw8/jfKIRI9cLD4cwZ/TxHDihSBHLm1G0rq+RJpPbt20eVKlXw8PCgcOHC+Pr6snDhwsQlUt9/DxUq6IXlVlYwa5YkUiLVkmQqlbC1tcXHx4cf7ZrR8dJiPFy/4O5dWUclxLuuVy+oXl0vKTIM8PbWGxEnh6tXr/Lll19SrVo1rl+/zsKFC6MSqwS5dk1/EICWLfXthu+9Z76AhbAQSaZSkYwZM9Lvkhdzvx3LSv+/qVChAjuXJ+NtO0KIJPfsmR6kuXdPt7t3h//8R1cHSC5Pnz5l2LBhlCxZkr/++ovBgwdz9uxZ2rVrh1VC7667fl0vLJ84UberVYO+fV/WixIiFZNkKhXqNdmTPXv24BJUjIqfN2Xm57J7jxDvigsXoEcPfWce6Jmwzz9P2mKbkZRS/PHHH5QsWZJffvmFZs2a4e/vz/Dhw8mSJcvbdxgWBpFrYz/8UG8IGLmnjRDvEEmmUikXFxfGbpjJiixufLd8BN26dSMkJMTSYQkhEmDSJBg6VD8vXVoX7+7VK3ljOHDgANWqVaNt27bkzZuX3bt3s2TJEmxtbRPeaf/+4OYGt2/r9rff6l2VhXjHSDKVin1UpQTtAtfi+WNfvObMZ3HOFhz8c7+lwxJCmCA4+OXzU6fg6NGXW9A5O+v1Ucnhxo0bdOzYkUqVKnHx4kXmzp3LgQMHqFGjRsI6PHgQIrfD+vpr+OMPkB0vxDtOkqlUztramlGjRuHVfyafhfjyn7Y/smXLFkuHJYSIw4YN8MEHcO6cbs+YoetFJedGB8HBwYwaNYoSJUqwZMkSBgwYwNmzZ+ncuXPC10Xduwc1aujpPICiRaF58+T9YEJYgCRT74i2Y7tybv0WDnx0m4YNG/Jr5+GEPQ+zdFhCiBeuX9e1KQHKlYPGjV+ug0qXLvniUEqxYsUK7O3t+fnnn6lfvz6nTp1izJgxZMuW7e07fPgQli7Vz3PlgpUrYZTUbRZpiyRT7xCXRpX43//+R4d6bflq/gTm5mtFYGCgpcMSIs0LDdULyX/4Qbfz54dFi/TATXI6cuQItWrVolWrVmTPnp3t27ezcuVKiiVmN+SJE6FtW/j3X91u1AgSkpQJkYpJMvWOyZIlC3P/9mKFswdjAw9RsWJFTpw4YemwhEhzzp59OduVLh3MmQPjx1smllu3btGtWzcqVKjAqVOnmDlzJocPH6Z27dpv35lS8Oef4Oen2999BwcO6B2WhUijJJl6B1lZW9H1yBQW7FpMUFAQmx374NVouKXDEiJN2bABRoyAy5d1u0kTXb08OT179oxx48ZRvHhxvLy86Nu3L+fOnaNHjx5YJ7TWwuPHeuuX6dN1O3t2PW8pRBomydQ7rEaNGuzbthdXqwfc2ridfv36ERoaaumwhHgn3bkDLVrA2rW63a2bTqQsUQlAKcXq1aspU6YM/fv3x83NjZMnTzJhwgRy5Mjx9h1euADDhulRKRsb2LULZs40e9xCpFaSTL3jipYphsNtH270cmLSpEm0Lv8Jp7bLtJ8Q5hJZqTxHDr1bSmQ7c2bImzf54zlx4gT16tWjRYsWZMiQgU2bNrF27VpKlCiR8E63bIFx417efli6tGwDI8QrDBVZ2CSZubi4qIORlXFFsli0cBGlOozmPRSPd8+meo3qlg5JiFStZ0/YsUPXibK21gM3lqoCcPfuXYYMGcJ///tfsmfPzvDhw+nRowfpEnKrYGgoTJkCpUrp+cmwMD30VqCA+QMXIpUwDOOQUsolpmMyMpWGtGvfjgfjxzAsVwFq16nN1ClTiQiPsHRYQqQaSsG6dRC52UCzZtC7N0S8+M/IEolUaGgokydPpnjx4syaNYtevXpx7tw5vv7664QlUqA/yPz5euEX6FEoSaSEiJWMTKVBDx8+pEOHDnyw9iENs2Sixpk/yPVhLkuHJUSKt2ePrkk5bx506mTpaGDDhg307duXM2fOUL9+fSZOnEiZMmUS1tk//+ipvEWLIEMGePAAcuY0b8BCpGIyMiWiyZEjB6tWraJi8ZJYP3lA7Qa18GhyrAAAH4ZJREFUOX/+vKXDEiLFUQp+/10nTwDVqsHq1Zbfq9ff3x93d3fc3d2JiIhg7dq1bNy4MeGJFMD9+7B/v15sDpJICfEWJJlKo6ysrOhydjbvrfuF6wHXqVu+NovaTbJ0WEKkCJED9oYBS5aAt/fLdrNmllt7/eDBAzw9PXFwcGDfvn1MmDCBEydO0KRJE4y3nWMMDdVzlFOn6naDBnD+vF5cLoR4K5JMpXGNGjfi0KFDDAhzoNXigQzp8hPh4eGWDksIi9m0CeztX96Vt3gxbN5s2ZjCwsKYMWMGxYsXZ9q0aXTp0oVz587Rr18/0qdP/3adRWaK6dLpquU3b+q2YejpPSHEW5NkSmBnZ0ebc4sY6tKREfNG07hxY25evWnpsIRINkFBL5OnQoX0WuvIds6clt2nd+vWrZQrV47evXvj6OjIkSNHmDlzJnny5Hn7zrZtA0fHlx9u9WoYPdq8AQuRBkkyJQB4/8P3+fWf/zJr1iyeb73Do8JV2Dx1raXDEiLJBQdD8eIwdKhuly6tyx0kpiyTOZw/f57mzZtTr149nj59ysqVK9m2bRuOjo5v31nYi03PCxSArFnh7l3dtpJfAUKYg/yXJKIYhkG3bt3o93N/rhh5+PJ7D+bPn2/psIQwu1u3Xq6DypQJfvkFPDwsGlKUwMBAfvjhB0qXLs327dsZM2YMJ0+epGXLlm+/LioiApo2hb59dbt0afD1hZIlzR+4EGmYJFPiDU2Gtcbp1nqca5aja+cuTC7SkaA7QZYOSwizmTQJOnbUSRXo4psuMd7wnHzCw8OZPXs2JUqUYMKECbRv355z584xYMAAMmbM+HadPX6sH62soGxZ+Ogj8wcshIgiyZSIUZ48edi4cSPDa3yL5+XfGejQgatXr1o6LCES5M7/t3ff4VFV6QPHvy81hpIAggWkiPSVFmQFC6CCgAoWkGYBEVBhRXAtK4qAoliwgAVRXIQfioAooiCrgFiIklAFpER6WZSSBIwkhLy/P87MzhACTEiZTOb9PM883HPvnTvn5M4Nb845971/wIABsHy5Kz/8MKxdC+edF9x6eS1ZsoRmzZrRv39/ateuTVxcHJMmTeL888/P/sHmz4cqVWDDBld+/nkYPDh3K2yMOYEFU+aUihUrxrDvXmXKQ2/zfymLaNq0KfM/nhfsahlzkiZN3CTxzK/Gjd32kiXh889hxQpXrlgx+HOiALZu3UrXrl1p3bo1Bw8eZPr06Xz33XfExMRk70Dp6b5J5c2aufwN55yT+xU2xmQpoGBKRNqLyEYRSRCRx7PYXlVEFovIShFZIyIdc7+qJljuevU+li1bRv1StWna/W7ebD6EYGXONyYrLVpA5gwBRYr45lmXLQtbt0K/fvlft6wcOXKEYcOGUa9ePebNm8eoUaPYsGED3bp1y/68KFVo3dqNW4KLFKdMgWrVcr3expisnTGYEpGiwJtAB6A+0ENEMmd1exKYoapNgO7AW7ldURNcdevWZcbCT/i2dDNei5tNly5dSE62eVSmYHjqqZNvTCtaFG67zXcjW3anHeWFjIwMPvjgA2rXrs1zzz1H165d2bhxI0899RTnZLcnyfvUAhHo29c938b+yDEmKALpmWoOJKjqFlVNA6YDnTPto0BZz3IUsCf3qmgKivNqns/tyfO4f+xg5syZw6sX9mDRG/ODXS0Tpg4fhpkzYdcud8d/mza+bSVKuF6o118PXrbyzJYuXcrll19O7969qVq1KrGxsUydOpUqVapk/2Dz57txym++ceU+fVzkGMyEWMaEsUCCqcqA/8zjXZ51/kYAd4jILmAe8I+sDiQi/UUkXkTi//jjj7Oorgk2EWHo0KHM/fccBvwZx8rBLzFz5sxgV8uEgaNHXTqDlStd+b//hdtvhy+/dOVx43xDfUWLut6qgmDnzp307NmTK664gt27dzN16tT/BVbZkpTkm1Tepg2MGhX8WxCNMUDuTUDvAUxW1SpAR2CqiJx0bFWdqKrNVLXZWWXvNQVGhztv4Pev/sPsmFRuv/12Hrv7nxw9cjTY1TKFiCpMneo6YcB1utxzjy8/1CWXwLJlbp233LevG+7r0wfO5ka43JSSksKIESOoU6cOn376KU899RSbNm3ijjvuoMjZJMts3x66dXM/mIgIePJJiI7O/YobY7JNzjSRWERaACNU9XpP+V8Aqvq83z7rgPaqutNT3gJcrqq/n+q4zZo10/j4+Jy3wARVWloaQwcNpf+7/2FHiQu4bMcMziso95ubkDNjhkuR5A2Q6tVzaZK8nZ+bNkHNmq7nKSt790L37vDxx8ELplSVjz76iMcee4xdu3bRrVs3XnjhBapld0K4qntQ4HXXubHKxYvdTPrs3ulnjMkVIrJcVbPsDg7kz6M4oJaI1BCRErgJ5p9n2mcHcK3nw+oBEYCN44WBEiVK8MbEN1je5nYmZKTStGlTYmNjg10tEyLmznWjVV4ffggTJvjKCxe6wMirdu1TB1Lg5k4tWRK8QCouLo4rrriCXr16UalSJb7//numT5+e/UAKYNEi6NDBRZjghvYskDKmQDpjMKWq6cAgYAHwK+6uvXUiMkpEOnl2exjoJyKrgY+A3mr3zoeVPoueZXT8O0RERDDxiseYFDOEjOMZwa6WKWAWLXJDcd7fDj/8AJMm+e64mzwZfvrJt/+FF4bG4+P27NnD3XffTfPmzdmyZQuTJk0iLi6OK6+8MnsH2rrVN6n8mmtcIHX77blfYWNMrjrjMF9esWG+wunggYP8WPk2IlOTmNyjAe+89w6RkZHBrpYJkvh4GDPG9Tadey68/z4MH+7mOl14oZtUXrJk6N6EdvToUV555RWee+45jh07xpAhQ3jiiScoW7bsmd+clbZt4bffYPPm03fBGWPyXU6H+YwJWPkK5emQ/DU/PNyZadOncU2zNsRO+z7Y1TL5JCEBbr7ZBVEAqalueds2V777bti50wVS4OZRh2IgparMmjWLevXqMWzYMNq1a8f69esZM2ZM9gKp1FQYPx68OdveeAO+/94CKWNCjAVTJtcVK1GMp19+mi+//JL7NkZwyR238OkHs4NdLZMHDh2CG2/0zWuKinLPvPM+QLhlSxdIee/gL1o0NIMnfytXrqR169Z07dqVsmXLsmjRImbPnk3NmjWzf7B16+DBB30z7OvUgcqZM88YYwo6C6ZMnunQoQP1po3m2QotubX3bYwYMYKMDJtHFcpU3WPfnn3WlaOi3CPh/vrLlStWdL1TN9zgyqEeOPnbt28f/fr1IyYmhvXr1zNhwgRWrFhBG/9soYH48UeYONEtN20Kq1e7iWTGmJBlwZTJU3/vfiXP75jOXXfdxVcjZ7LgnNZsW7k12NUyZ+A/lbJPH7jjDrcs4lIbeafBFSkCsbHQu3e+VzHfpKam8tJLL1GrVi0mT57MkCFD2Lx5MwMGDKDo2QzHvfMOvPACpKW5csOGuVthY0y+s2DK5LnIyEgmT57MPVfeRLW03dzU6SZWrVoV7GoZP8eO+ZYff/zExNo1asDFF/vKU6bA0KH5V7dgUVXmzJlDgwYNePTRR2nVqhXr1q1j7NixRGcnWebBgzB4sG/i2CuvwJo1Jz+Z2RgTsiyYMvlCROj//RgOLZ7CIU2kxd9b8G7X58/8RpMn/vzT1/v0+utueC411ZUbNIDWrcE7Ijt8+Im5oMLB2rVradu2LTfffDMlS5ZkwYIFzJ07l9q1a2f/YCkp8O9/u6Sb4G5rLFUqdytsjAkqC6ZMvrqi9RWsWLGCIeVvpN+sJxh21SDSvMMdJs8kJ/uCpU8+cUN1W7a4cuPGMGCAb97TnXfC2LGhkd8pt+3fv58HHniARo0asWLFCsaPH8+qVato165d9g40Y4av+65KFdi+3Y2XGmMKpTD8dWmCrVKlSgzfPJUXmw7iuR/epFWrVuzYtiPY1SpUDh2CxES3vGwZlCvnywXZuDE88ojL7wTQqpWbwhPOj3k7duwYr732GrVq1WLixIkMHDiQhIQEBg0aRPHixbN/wHXr4LvvXK8UuBNgjCm0LJgyQRFROoJHl49n5syZ7F+1j6Qa1zFryPvBrlbIOnAAdu1yy7//DhUquASZ4J5tN2wY1KrlyjVrwnPPuQ4TA/Pnz+fSSy9lyJAhNG/enDVr1jBu3DjKly8f+EH27IGbboKlS135iSdcFGsJa40JCxZMmaDq0qUL742bSKqU4OnXn+fll1/GnkR0Zvv2wfr1bjkjAy65BJ55xpUrVXJznNu2deXISDfn6Wym+xRmGzZsoGPHjnTs2JGMjAzmzp3LV199Rf369QM/iPe7GhXlMpdv3+7KJUuG5zipMWHKrnYTdK36XUftQ0upd2sjHnnkEYbXupu9m/YEu1oFyp49Lj2R1w03wAMPuOUiRdzjWgYM8G1/6CG49NL8rWOoOHToEA899BCXXnopS5cuZezYsaxdu5Ybb7wRyU5irLfegnbtXEBVqhT88gv06JF3FTfGFFgWTJkCoWxUWWbOnMkL/Z/l6d8+ZGrje9mwYUOwqxU0e/bA55/7yo88Al26+DpCXnwRXn7Zt71bN5f/0Zxaeno6b731FrVq1WL8+PH07duXzZs3M3ToUEoEmqbg+HHfbY7nnANly8Lhw65sj4AxJmxZMGUKDBHh0XeG8dnQd3k9cjWXXXYZH3/wcbCrlS/++1+YPBnS0135vffcM+4OHXLlxx6DL77w7X/NNSfmgjKn980339CkSRMGDhxIw4YNWblyJRMmTKBixYqBH2TPHmjSxPfsnN693a2RZ/tQY2NMoWHBlClwuoztw8+rfqZxnUbU7D2S96reRbo3yigk9u+Hd991QRTAkiXuzvmVK13ZuxwV5coNG0JMTOF6PEt+SEhIoHPnzrRt25aUlBRmz57NwoULaZidrON//un+Pf98N4vfe1LsZBhjPCyYMgVSlSpVmLdgPpvL1WHOzg20a9eO33//PdjVOmuHD8Pbb7vE1+DuvOvfH7791pU7dHAPCI6JceWLLoJGjWwOc3ZMmzaN6tWrU6RIEapWrcoNN9xA/fr1WbRoEWPGjGHdunXccsst2ZsXNWaMux3yr7/cyfjkE+jYMe8aYYwJSfar2hRYZSqUocfBT+kyeSCxsbH8s2YvPh/2YbCrFZDjx2H8eFiwwLdu0CD48ku33LAh/Pqrm+sEbqSoQQMLns7WtGnT6N+/P9u3b0dV2blzJ/PmzaNly5Zs3ryZxx57jIiIiMAOdviwL4NpixZw222+8VdjjMmCBOs29GbNmml8fHxQPtuEnvif4znn8p7soQwJb97Lffffl70ehnwwYQIULw59+7qJ4tWqubvu3n7bbd+1CypXttGhU0lPTyc5OZnExMSTXklJSVmu927bsWNHlik1qlWrxjbvM/EC8ccf7jbIBx90uaKMMcZDRJarapazVYvld2WMORvN/t6MrfFfMWHww8we+ABxS5bx8itjKV85G4kVc9nUqZCQACNHuvInn7j0Qn37uoBp9eoTE18X9iSZx44dO2PQc6ptiYmJHDly5IyfERUVRXR09P/+rV69OtHR0UyZMiXL/XfsCCCzvqrLEXXJJe4hhQMG+JJ0GWNMAKxnyoSU48ePM2rUKGqP+oa6RZIos2Y2tRvkTzbKzz6D2bPB+//2oEHw008QF+eCp7/+cnfLh6rU1NSAe4KyWp/ifXTKKRQpUuR/QZD3lbl8um1lypSh6CnSD1SvXp3t3oSZfgLqmRoxwuWZSEhwk8yNMSYL1jNlCo2iRYsycuRI/m9TFNNnzWPSVZfz4Ycf0r59+1z/rMWL4fnnXY9TmTJumG75ckhKcjd0vfYaFPO7goIZSKkqR48ezfbwmH/56NGjp/2MokWLnhTwXHDBBQEHRaVLl6ZIHk0KGz16NP379z8hoIuMjGT06NFZv2H7dnfCKlWCXr1cj1SFCnlSN2NM4Wc9UyZk/fbbb9x6660UX3OMB2pczl0bJlKsxNn/fbB6tcscPm6cmzbzzTfw8MMwYwbUqeNGg/JqvpOqkpKSclZDZN5taWlpp/2M4sWLU65cuWz1Bvm/IiMjC9w8NX/Tpk1j2LBh7Nixg6pVqzJ69Gh69ep18o7JyW7MtWdPN9HNGGMCcLqeKQumTEhLSUlhetU7ufbAzwxt05xJs98nOjo6oPfu2ePuphs6FG65xXVW3Hqr63G66qrs1UNVOXz48FnPGUpKSjpjLq2IiIizHiKLjo4mIiKiQAdDeSotzSXz8s6F+ugjuPJKl4PCGGMCYMN8ptCKjIyk976ZvPH0OGaPbs3sclH8nVhe4HG68TH7OJ/GjV0CzGPHoHVr6NTJZRSvWNGlIhCBjIwMoqKS+eQTF9wsWZK9oCgpKYkM72NGTlNX/+DmvPPOo06dOgEFRVFRUYHf2m9O9tJL8NRTsHGjS7xpz9AzxuQi65kyhcZll8XRMv57XuZRipLB29zPG9xP2ch9tP/ndyQmJrLms5aUKrWWwxW/IykpiRL795OYnEzCkSOoKuWBdCDZc8wIIANIA8qUKUNU2bJElyuX7SGyqKgoihcvHqwfTfjYuxe6d3ePfNm2zSXwql/fPZcnNtZlRw3X3jljTI6crmcKVQ3KKyYmRo3JTZUrN9NL+UnTKKYK+ifn6EJa6U9cqIBGRUXp0pIlNb5UKb366qu1U6dOuqliRd1w0UU6fPhwfeWVV/T36tV1b9OmunDhQl2+fLkerVtXUzt21GPHjrkPqVNHtVs334c2aKDar5+v3KiR6pAhvnJMjOrw4b5yixaqL77oK199teqbb/rK112nOnmyW05PV+3USXXGDFc+elS1Z0/VuXNd+cgR99kLF7pycrLqQw+p/vijKx86pPrkk6rLl7vygQOqzz+vunatK+/frzpunOrmzb7y+++rbt/uK3/8serevb7yl1+6f1VVDx5U/fZb1cRE3+fFxbl6qaomJan++qurt7e+O3aopqX52nPggGunt72pqaoZGXrW7r9ftUgR1f79VStUOPFcGWNMDgDxeoqYxvItm0Jjz57l3McjHPd8rYtwnP0kMZQ9pKenk5iYSIsFC4j54guWLFnCnDlzqPXBB9R57z1GjhzJkCFDqPjSS5w/ciTXXHMNTZs2peSjj1Li3nsp5r1tb+BAlxHbq2dPaNPGV27f3j0M16t5c7j4Yl+5Rg0491xfuXRpl5wK3Az3v/5y45EAGRluIldioiunp8OyZa73Bdw8oC++gK1bXfnPP2HSJJdaHdz7Ro/2PcNm/37417/cTHv3A3PJKb3lbdvgnntg1SpX3rjRTSrzbv/lF5eF1Hu8uDg3brpunSt//z1cdpnv87/+GurVc8cB+PxzqFoVtmxx5enT3R103lxQkya5n4W3fW+9BSVKuHoDvPGGS9x1+LArjx/vJpKnprrys8+6CeUZGS4J2AcfuCdGG2NMHgtomE9E2gOvA0WB91R1TKbtrwLe/1EigUqqetpZwDbMZ3LbZVWq8N3u3fhnKEgBWlWpQtzOncGqVsGRkeECj+LFXU6H9HQXcJUuDRERbtvevS5dQGQkpKS4QK1aNbdPcjJs2AB167rhs4MHXeAVE+NyRezb54K9q66C6GjYuRN++MENrUVHu8SYixdD165u//XrXcDVp4873sqVMG8eDB7sPi82FubOhSefdPVZtAjmzIEXX3RB17x5Lm/FO++49lx/vbsFMyPDBWH33gtvvhnsn7oxppDI0d18IlIU2AS0BXYBcUAPVV1/iv3/ATRR1XtOd1wLpkxu23TddVRduBD/adpHge3XXkudb74JVrVMfti71/UA+ufKOucc1wtmiTiNMbngdMFUIMN8zYEEVd2iqmnAdKDzafbvAXyU/WoakzO1Dxwg8/1uEUCdAweCUR2Tn555xvVI+Tt+3K03xpg8FkhqhMqA/xjJLuDvWe0oItWAGsCiU2zvD/QHqFq1arYqaswZrVwZ7BqYYImNdXPI/KWlwdKlwamPMSas5Haeqe7ALFU9ntVGVZ0ITAQ3zJfLn22MCVcWSBtjgiiQYb7dgH+a4CqedVnpjg3xGWOMMSaMBBJMxQG1RKSGiJTABUyfZ95JROoC5YDY3K2iMcYYY0zBdcZgSlXTgUHAAuBXYIaqrhORUSLSyW/X7sB0DSTXgjHGGGNMIRHQnClVnQfMy7RueKbyiNyrljHGGGNMaLAM6MYYY4wxORC0Bx2LyB/A9nz4qHOB/fnwOQWRtT18hXP7w7ntEN7tt7aHr/xofzVVrZjVhqAFU/lFROJPlbG0sLO2h2fbIbzbH85th/Buv7U9PNsOwW+/DfMZY4wxxuSABVPGGGOMMTkQDsHUxGBXIIis7eErnNsfzm2H8G6/tT18BbX9hX7OlDHGGGNMXgqHniljjDHGmDxjwZQxxhhjTA6EbDAlIu1FZKOIJIjI41lsLykiH3u2/ywi1f22/cuzfqOIXJ+f9c4NAbR9qIisF5E1IrJQRKr5bTsuIqs8r5OesRgKAmh/bxH5w6+d9/ptu1tENnted+dvzXMugLa/6tfuTSKS6LctpM+9iLwvIr+LyNpTbBcRGef52awRkaZ+20L6vENA7e/lafcvIrJURBr5bdvmWb9KROLzr9a5I4C2txaRJL/v93C/bae9Zgq6ANr+iF+713qu8/KebSF93gFE5CIRWez5P22diAzOYp/gX/uqGnIvoCjwG3AxUAJYDdTPtM8DwATPcnfgY89yfc/+JYEanuMUDXabcrntbYBIz/L93rZ7ykeC3YZ8aH9v4I0s3lse2OL5t5xnuVyw25Sbbc+0/z+A9wvRub8aaAqsPcX2jsB8QIDLgZ8Lw3nPRvtbetsFdPC231PeBpwb7DbkYdtbA19ksT5b10xBfJ2p7Zn2vQlYVFjOu6cNFwBNPctlgE1Z/M4P+rUfqj1TzYEEVd2iqmnAdKBzpn06Ax94lmcB14qIeNZPV9VUVd0KJHiOFyrO2HZVXayqKZ7iT0CVfK5jXgrk3J/K9cDXqnpQVQ8BXwPt86ieeSG7be8BfJQvNcsHqvodcPA0u3QGpqjzExAtIhcQ+ucdOHP7VXWpp31QyK77AM79qeTk90WBkM22F6prHkBV96rqCs/yYeBXoHKm3YJ+7YdqMFUZ2OlX3sXJP9z/7aOq6UASUCHA9xZk2a1/X1zE7hUhIvEi8pOI3JwXFcxjgbb/Nk937ywRuSib7y2oAq6/Z2i3BrDIb3Won/szOdXPJ9TP+9nIfN0r8B8RWS4i/YNUp7zWQkRWi8h8EWngWRc2515EInGBwid+qwvVeRc3XacJ8HOmTUG/9ovlxUFNwSAidwDNgFZ+q6up6m4RuRhYJCK/qOpvwalhnpkLfKSqqSIyANdDeU2Q65TfugOzVPW437pwOPdhT0Ta4IKpK/1WX+k595WAr0Vkg6fHo7BYgft+HxGRjsBnQK0g1ym/3QT8qKr+vViF5ryLSGlcoPiQqiYHuz6ZhWrP1G7gIr9yFc+6LPcRkWJAFHAgwPcWZAHVX0SuA4YBnVQ11bteVXd7/t0CfIuL8kPJGduvqgf82vweEBPoewu47NS/O5m6+wvBuT+TU/18Qv28B0xEGuK+851V9YB3vd+5/x34lNCa2nBGqpqsqkc8y/OA4iJyLmF07jn9NR/S511EiuMCqWmqOjuLXYJ/7QdjQllOX7getS24YQzvpMIGmfYZyIkT0Gd4lhtw4gT0LYTWBPRA2t4EN+myVqb15YCSnuVzgc2E3mTMQNp/gd/yLcBPnuXywFbPz6GcZ7l8sNuUm2337FcXN/FUCtO599S9OqeehHwDJ05CXVYYzns22l8VNwe0Zab1pYAyfstLgfbBbksut/187/cdFzDs8HwPArpmCvrrdG33bI/CzasqVQjPuwBTgNdOs0/Qr/2QHOZT1XQRGQQswN2t8b6qrhORUUC8qn4OTAKmikgC7kvW3fPedSIyA1gPpAMD9cShkAItwLa/BJQGZro59+xQ1U5APeAdEcnA9UqOUdX1QWnIWQqw/Q+KSCfc+T2Iu7sPVT0oIs8AcZ7DjdITu8QLtADbDu67Pl09v008Qv7ci8hHuLu2zhWRXcDTQHEAVZ0AzMPd1ZMApAB9PNtC+rx7BdD+4bh5oW95rvt0VW0GnAd86llXDPhQVb/K9wbkQABt7wLcLyLpwF9Ad8/3P8trJghNOGsBtB3cH43/UdU//d4a8ufd4wrgTuAXEVnlWfcE7o+HAnPt2+NkjDHGGGNyIFTnTBljjDHGFAgWTBljjDHG5IAFU8YYY4wxOWDBlDHGGGNMDlgwZYwxxhiTAxZMGWMKPBGJFpEHPMsXisisYNfJGGO8LDWCMabA8zyT6wtV/VuQq2KMMScJyaSdxpiwMwao6Unatxmop6p/E5HewM24DM+1gJdxma7vBFKBjp7EfTWBN4GKuKR+/VR1Q/43wxhTGNkwnzEmFDwO/KaqjYFHMm37G3ArcBkwGkhR1SZALHCXZ5+JwD9UNQb4J/BWvtTaGBMWrGfKGBPqFqvqYeCwiCQBcz3rfwEaep423xLf45XAPZvTGGNyhQVTxphQl+q3nOFXzsD9jisCJHp6tYwxJtfZMJ8xJhQcBsqczRtVNRnYKiJdAcRplJuVM8aENwumjDEFnqoeAH4UkbXAS2dxiF5AXxFZDawDOudm/Ywx4c1SIxhjjDHG5ID1TBljjDHG5IAFU8YYY4wxOWDBlDHGGGNMDlgwZYwxxhiTAxZMGWOMMcbkgAVTxhhjjDE5YMGUMcYYY0wO/D9M1hfLVPU4uQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" } } ] }, { "cell_type": "markdown", "metadata": { "id": "WIlW946pPMxy" }, "source": [ "## Local Error \n", "The error for the 1 step Adams Moulton is:\n", "\\begin{equation}\\frac{y_{n+1}-y_{n}}{h}=\\frac{1}{2}[f(t_{n+1},w_{n+1})+f(t_{n-1},w_{n-1})] +\\frac{-h^2}{12}y^{(3)}(\\eta),\\end{equation}\n", "where $\\eta \\in [t_{n-1},t_{n+1}]$.\n", "The error for the 2 step Adams Bashforth is:\n", "\\begin{equation}\\frac{y_{n+1}-y_{n}}{h}=\\frac{1}{2}[3f(t_{n},w_{n})-f(t_{n},w_{n})] +\\frac{5h^2}{12}y^{(3)}(\\xi),\\end{equation}\n", "where $\\xi \\in [t_{n-1},t_{n+1}]$.\n", "\n", "Rearranging the equations gives \n", "\\begin{equation}\\frac{w_{n+1}-w_{n+1}^0}{h}= h^2[\\frac{y^{(3)}(\\eta)}{12}+\\frac{5y^{(3)}(\\xi)}{12}]\\approx h^2[\\frac{6}{12}y^{(3)}(\\xi)].\\end{equation}\n", "Making the assumption that \n", "\\begin{equation}y^{(3)}(\\xi)\\approx \\frac{12}{6h^2}\\frac{w_{n+1}-w_{n+1}^0}{h},\\end{equation}\n", "\\begin{equation}y^{(3)}(\\xi)\\approx \\frac{12}{6h^2}\\frac{w_{n+1}-w_{n+1}^0}{h}.\\end{equation}" ] }, { "cell_type": "code", "metadata": { "id": "hC10-rujPMxy", "outputId": "31399ea8-bd9b-442c-934b-ccd595046d1e", "colab": { "base_uri": "https://localhost:8080/", "height": 0 } }, "source": [ "\n", "d = {'time t_i': t, 'Adams Predictor w0': w0,\n", " 'Adams Corrector':w,'Exact (y)':y,'|w-y|':np.round(np.abs(y-w),5),'|w0-w|':np.round(np.abs(w0-w),5),'estimate LTE':np.round(1/(6*h)*abs(w0-w),5)}\n", "df = pd.DataFrame(data=d)\n", "df" ], "execution_count": 11, "outputs": [ { "output_type": "execute_result", "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
time t_iAdams Predictor w0Adams CorrectorExact (y)|w-y||w0-w|estimate LTE
00.01.0000001.0000001.0000000.000000.000000.00000
10.50.7130610.7130610.7130610.000000.000000.00000
21.00.8032650.7089800.7357590.026780.094290.03143
31.50.9805100.9116070.9462600.034650.068900.02297
42.01.2801471.2386691.2706710.032000.041480.01383
\n", "
" ], "text/plain": [ " time t_i Adams Predictor w0 ... |w0-w| estimate LTE\n", "0 0.0 1.000000 ... 0.00000 0.00000\n", "1 0.5 0.713061 ... 0.00000 0.00000\n", "2 1.0 0.803265 ... 0.09429 0.03143\n", "3 1.5 0.980510 ... 0.06890 0.02297\n", "4 2.0 1.280147 ... 0.04148 0.01383\n", "\n", "[5 rows x 7 columns]" ] }, "metadata": {}, "execution_count": 11 } ] }, { "cell_type": "code", "metadata": { "id": "rDTUnci_PMxz" }, "source": [ "" ], "execution_count": 11, "outputs": [] } ] }